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ABSTRACT 

This  report  describes  the  development  of  a  computer  code  to  solve  the  two-dimensional 
boundary  layer  equations  in  direct,  inverse  or  mixed  direct/inverse  mode  for  airfoil  flows. 
The  solution  algorithm  uses  body  conformal  coordinates  and  a  relaxation  scheme  with 
flow-dependent  operators.  The  code  incorporates  two  methods  to  sweep  the  flow  field. 

Topics  include  a  description  of  the  code  structure,  its  input  requirements,  boundary 
condition  and  flow  field  initialization  and  various  software  tools  required  by  the  main 
algorithm.  Finally,  verification  results  are  presented. 

This  work,  done  in  cooperation  with  the  University  of  Toronto,  seeks  the  development  of 
a  boundary  layer  code  compatible  with  the  NASA  Ames  ARC2D  Navier-Stokes  code. 
The  two  codes  will  be  used  in  a  study  of  the  Fortified  Navier-Stokes  concept. 


RESUME 

Ce  rapport  ddcrit  le  ddveloppement  d’un  logiciel  pour  rdsoudre  les  dquations  des  couches 
limites  appliquees  a  l’dcoulement  autour  d’un  profil  d’aile.  Les  Equations  sont  exprimees 
en  coordonndes  conformes  au  corps  dtudid;  une  mdthode  de  relaxation  avec  des 
opdrateurs  s’ajustant  k  l’dcoulement  local  ainsi  que  des  mdthodes  de  mode  direct,  inverse 
et  mixte  direct/inverse  sont  employees.  Enfin,  deux  mdthodes  de  balayage  du  champ 
sont  utilisdes. 


Les  sujets  suivants  sont  abordes:  la  structure  du  logiciel,  les  entities  necessaires  a  son 
operation,  les  conditions  limites  et  l’initiation  des  variables.  De  plus,  differents  outils 
numeriques  requis  par  le  logiciel  principal  sont  ddcrits  et  des  rdsultats  ddmontrant  la 
validitd  du  lociciel  sont  presents. 


Ce  travail,  mend  en  collaboration  avec  l’Universite  de  Toronto,  vise  k  dlaborer  un 
programme  de  couche  limite  compatible  avec  le  programme  Navier-Stokes  ARC2D  mis 
au  point  par  le  centre  de  recherche  Ames  de  la  NASA.  Les  deux  programmes  serviront  k 
dtudier  le  concept  fortifid  de  Navier-Stokes. 
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INTRODUCTION 


Background 

In  the  mid  1980’s,  researchers  at  NASA  Ames  proposed  the  Fortified  Navier-Stokes 
(FNS)  concept,  a  procedure  by  which  approximate  solutions,  or  any  other  known 
information,  is  used  to  improve  the  accuracy  and  efficiency  of  a  Navier-Stokes  code. 
Between  1985  and  1988,  Steger  and  Van  Dalsem'1"2"3'  explored  the  FNS  concept  by 
combining  a  Navier-Stokes  code  with  a  boundary  layer  code.  The  boundary  layer 
equations  are  solved  in  a  body  conformal  frame  of  reference  while  the  Navie-  Stokes 
solver  uses  Cartesian  velocity  components  with  generalized  curvilinear  coordinates. 
Since  the  two  procedures  use  different  grids  and  velocity  components,  going  from  one 
solution  to  the  other  involves  interpolation  and  rotation  of  the  velocity  data.  Despite 
these  differences,  Steger  and  van  Dalsem  still  improved  the  convergence  of  their 
Navier-Stokes  code  by  an  order  of  magnitude.  Their  work  laid  a  solid  foundation  for 
further  investigations  of  the  FNS  concept. 

The  purpose  of  the  current  project  is  to  investigate  the  range  of  validity  of  the  FNS 
approach  as  a  function  of  factors  like  Mach  number,  Reynolds  number,  and  angle  of 
attack.  Three  phases  are  planned: 

i.  the  development  of  a  boundary  layer  solver  compatible  with  the  NASA  Ames  code 
ARC2D, 

ii.  integration  of  the  boundary  layer  and  Navier-Stokes  codes,  and 

iii.  parametric  studies  using  the  FNS  code. 

This  report  presents  partial  results  for  phase  i)  above.  It  describes  the  development  of 
GBL,  a  boundary  layer  research  code  to  study  the  Steger  and  Van  Dalsem  algorithms. 
GBL  has  the  flexibility  to  suppoit  changes  to  the  formulation  of  the  equations,  the 
boundary  conditions,  or  to  the  various  iterative  schemes  used  to  solve  the  equations. 

At  this  point,  it  is  useful  to  review  the  boundary  layer  work  of  Steger  and  Van  Dalsem. 


Previous  Generalized  Boundary  Layer  Work 

Several  boundary  layer  methods  (e.g.  Cebeci'4')  employ  a  coordinate  transformation,  or 
"stretching"  of  the  governing  equations  prior  to  formulation  of  the  finite  difference 
equations.  The  purpose  is  to  scale  the  growth  of  the  viscous  layer  so  that  the  new 
equations,  expressed  in  terms  of  similarity  variables,  may  be  solved  on  an  approximately 
rectangular  grid.  Van  Dalsem  and  Steger5'  took  a  different  approach.  They  developed  a 
scheme  to  solve  the  two-dimensional,  steady  state  boundary  layer  equations  in  body- 
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conformal  coordinates  (s,n),  which  amounts  to  solving  the  equations  in  physical  -  instead 
of  similarity  -  variables.  An  adaptive  grid  must  then  be  used  to  concentrate  a  sufficient 
number  of  computational  points  within  the  boundary  layer.  Solution  of  the  equations  on 
a  grid  with  uneven  spacing  results  in  complicated  finite  difference  equations,  but  this  is 
avoided  by  transforming  the  equations  to  a  computational  domain  with  uniform  grid 
spacing  so  that  standard,  unweighted  differences  can  be  used. 

The  scheme  of  Van  Dalsem  and  Steger  is  noteworthy  for  two  reasons.  First,  it  assumes 
the  equations  are  only  weakly  coupled  when  the  pressure  is  given.  Each  equation  can 
then  be  solved  independently  in  a  sequential  scheme.  They  also  chose  their  finite 
difference  operators  to  obtain  bidiagonal  or  tridiagonal  systems  of  equations  which  are 
solved  efficiently  with  the  Thomas  algorithm. 

Second,  their  scheme  avoids  use  of  complex  space  marching  procedures  within  separated 
flow  regions.  Instead,  the  scheme  marches  in  a  single  direction  while  flow-dependent 
finite  difference  operators  adapt  locally  to  respect  the  parabolic  nature  of  the  boundary 
layer  equations. 

The  earlier  work  of  Steger  and  Van  Dalsem151  used  a  predictor-corrector  finite  difference 
approximation  to  the  two-dimensional,  steady-state  equations.  A  mix  of  upwind  and 
central  operators  is  used  to  differentiate  the  convective  terms  within  the  momentum  and 
energy  equations.  The  choice  of  which  operator  to  use  depends  on  the  local  value  of 
u/ue.  This  algorithm  was  later  expanded161  to  the  three-dimensional  unsteady  form.  The 
formulation  was  also  simplified  to  a  single-step  process  and  only  the  upwind  operator, 
dependent  on  the  sign  of  the  contravariant  velocity  U  or  W,  was  retained  for  the 
convective  terms.  The  reader  should  be  aware  that  Steger  and  Van  Dalsem  do  not  use 
the  time  variable  in  a  true  time  fashion,  since  the  method  still  assumes  that  the  pressure 
distribution  is  fixed  and  given.  Instead,  the  "time-like"  variable  is  used  to  relax  the 
equations. 

The  above  review  would  not  be  complete  without  mention  that,  in  1988,  Steger  and  Van 
Dalsem^71  developed  a  very  different  approach  to  solving  the  boundary  layer  equations. 
Having  recognized  the  lack  of  a  boundary  layer  procedure  directly  compatible  with  their 
Navier-Stokes  algorithm,  they  recast  the  boundary  layer  equations  in  a  form  using  the 
Cartesian  velocity  components  instead  of  the  usual  body  conformal  components.  The 
new  equations  can  be  solved  on  a  Cartesian  grid,  thus,  increasing  commonality  with 
Navier-Stokes  codes. 


GBL  Program  Structure 

GBL  uses  the  1986  body  conformal  version  of  the  Steger  and  Van  Dalsem  boundary 
layer  equations  with  time-like  relaxation.  It  was  selected  over  the  1988  version  because 
it  is  simpler  to  implement  and  provides  a  framework  from  which  to  study  the  behavior  of 
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the  algorithm.  Scaling  of  the  equations,  however,  was  modified18^  to  increase  their 
compatibility  with  Pulliam’s  ARC2D^  Navier-Stokes  code.  The  latter  will  be  used  as 
the  Navier-Stokes  components  of  the  FNS  code. 

The  infrastructure  of  GBL  is  designed  in  four  general  parts:  input,  initialization, 
computations  and  output.  In  addition,  specialized  graphics  tools  were  developed  to 
control  program  execution  and  to  examine  data  during  the  development  cycle.  However, 
these  tools  are  machine  specific  and  their  description  is  outside  the  scope  of  this  work. 
For  these  reasons,  the  graphics  tools  are  not  described  in  this  report. 

Program  inputs  are  described  in  the  PROGRAM  INPUTS  section.  They  include 
parameters  specifying  the  freestream  conditions  as  well  as  switches  to  control  program 
execution.  In  addition,  the  user  may  elect  to  run  one  of  the  built-in  test  cases  to  verify 
the  implementation  of  a  new  algorithm,  or  to  run  an  airfoil  flow  field  calculation.  The 
user  may  also  direct  GBL  to  generate  its  own  grids  for  the  test  cases  or  select  one  of  two 
grid  generation  schemes  for  the  airfoil  case  (an  "adaptive"  grid  or  a  subset  of  the  Navier- 
Stokes  grid).  Various  tools  are  required  by  the  input,  initialization  and  computation 
routines  as  well  as  for  output  processing.  These  are  described  in  the 
COMPUTATIONAL  TOOLS  section. 

Four  initialization  procedures  may  be  used  by  GBL:  two  use  ARC2D  data  to  initialize  the 
flow  about  an  airfoil  while  the  others  are  test  cases  for  flow  over  a  flat  plate.  A 
description  of  these  procedures  is  provided  in  the  INITIALIZATION  OF  VARIABLES 
section.  In  the  SOLUTION  SCHEME  section,  the  algorithm  and  overall  program  flow 
along  with  the  application  of  boundary  conditions,  sweeping  schemes  and  determination 
of  the  inverse  mode  region  are  described. 

The  CODE  TESTING  section  present  code  validation  results  while  conclusions  are 
presented  in  the  last  section. 
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PROGRAM  INPUTS 

The  user  selects  the  run  parameters  for  GBL  by  editing  the  control  file  inplrun.sw  which 
contains  sufficient  information  for  GBL  to  proceed  automatically,  whether  it  is  a  test  case 
or  airfoil  computations.  The  parameters  are  divided  into  two  general  categories:  essential 
parameters,  which  are  always  required  by  the  program,  and  test  parameters  which  are 
used  to  verify  new  algorithms,  as  listed  below. 

Essential  parameters: 

[1]  Freestream  parameters,  used  to  compute  all  freestream  variables  required  to  scale 
the  equations. 

[2]  Turbulence  model  switch  and  transition  points  for  lower  and  upper  surfaces  of  the 
airfoil. 

[3]  Grid  parameters  which  include  a  switch  to  chose  between  a  flat  plate  grid 
generated  internally  or  an  airfoil  grid  read  from  ARC2D.  The  flat  plate  grid 
generator  requires  user  supplied  specifications  while  for  airfoil  cases,  the 
boundary  layer  grid  is  generated  automatically  from  the  Navier-Stokes  grid.  Grid 
generation  is  discussed  in  further  details  later  in  this  section. 

[4]  Time  step  to  use  for  the  time-like  relaxation  scheme  and  total  number  of  cycles 
for  the  run. 

[5]  Convergence  criteria  to  use  for  computation  of  the  residuals. 

Test  parameters: 

[1]  Mode  (direct/inverse)  to  be  tested. 

[2]  Direction  for  sweeps  of  the  computational  field  during  the  solution  process. 

[3]  Switches  to  control  execution  of  individual  equations  during  the  solution  process. 

[4]  Self-similarity  parameter  for  test  cases  involving  Falkner-Skan  flow. 

[5]  Control  over  the  %  range  of  the  computation;  particularly  useful  when 
investigating  behavior  of  the  algorithm  over  a  specified  flow  region,  e.g.  about  the 
airfoil  nose  or  in  its  wake. 

[6]  Control  over  the  extent  of  the  region  where  the  inverse  formulation  is  used  as 
well  as  specification  of  parameters  used  to  update  the  inverse  forcing  function  (to 
be  explained  later). 

[7J  Control  over  individual  boundary  conditions  for  the  momentum  and  energy 
equations. 
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[8]  Selection  of  an  additional  relaxation  factor  (extra  to  the  time-like  terms)  which 
may  be  applied  to  the  solution  of  individual  equations. 

A  sample  control  file  is  given  in  appendix  A  along  with  a  brief  description  of  each 
parameter. 

GBL  supports  flow  field  initialization  for  several  test  and  airfoil  cases.  The  test  cases 
involve  Falkner-Skan  or  Klineberg  flow.  These  flows  are  stricly  laminar,  however,  the 
Falkner-Skan  case  may  also  be  used  with  turbulence  to  simulate  turbulent  flat  plate  flow. 
For  either  case,  GBL  generates  a  flat  plate  boundary  layer  grid  according  to  user 
specifications.  Airfoil  cases  require  the  user  to  provide  a  Navier-Stokes  grid  and  the 
corresponding  state  vector  from  ARC2D.  GBL  then  geneiates  a  compatible  boundary 
layer  grid  and  initializes  flow  field  variables  from  the  ARC2D  data. 

The  remainder  of  this  section  provides  details  about  computation  of  the  freestream 
parameters,  generation  of  the  flat  plate  and  airfoil  boundary  layer  grids  and  their 
transformation  to  body  conformal  coordinates. 


Freestream  Parameters 

The  freestream  parameters  are  used  to  scale  flow  field  variables.  The  user  provides  four 
parameters:  (i)  the  Reynolds  number  ( Re «,),  (ii)  the  Mach  number  (M^),  (iii)  the 
freestream  temperature  ( T «,)  and  (iv)  a  reference  length  (/)•  All  other  freestream 
parameters  are  computed  from  these  four  values  with  the  assumption  that  the  fluid 
behaves  like  a  perfect  gas  and  that  the  flow  is  isentropic.  The  computations  proceed  in 
the  following  order. 

The  freestream  speed  of  sound  is  calculated  from  the  freestream  temperature,  the  ratio  of 
specific  heat  at  constant  pressure  to  specific  heat  at  constant  volume,  y  (1.4  for  air),  and 
the  gas  constant  R  (287.3  m2/K-sec 2  for  air). 


a„  =  (yRTjA 


(2-1) 


Having  determined  a„,  the  freestream  velocity  u „  is  obtained  from  the  Mach  number 
definition 


=M00aeo 


(2-2) 


The  freestream  temperature  is  also  used  with  Sutherland’s  approximation*10*  (p.  312)  and 
reference  values  of  the  temperature  and  viscosity  at  sea  level,  as  defined  in  the  standard 
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atmosphere  (Tsl  =  288. 16K,  [isi  =  1.78935x10  5  N-sec/m2),  to  compute  freestream 
viscosity 


= 


1.458xlQ~6  Tt2 
T  S\ 


(2-3) 


where  Sj  =  1 10.4  K.  Freestream  density  is  obtained  from  the  definition  of  the  freestream 
Reynolds  number  and  the  reference  length 


P-  = 


Reo4U 

u^l 


(2-4) 


while  the  kinematic  viscosity  is  simply  the  ratio  of  viscosity  to  density 

Hoo 

Vo,  = - 

Poo 


(2-5) 


and  freestream  pressure  is  computed  from  the  perfect  gas  relation 

Poo  =  Poo  RTm  (2-6) 


Finally,  freestream  total  enthalpy  is  computed  by  substituting  freestream  values  in  its 
definition 


H„  =  cpT  +  (2-7) 

Note  that  the  normal  velocity  x>„  is  omitted  from  the  above  definition  of  total  enthalpy 
because  of  the  body  conformal  assumption  used  to  derive  the  equations. 


Flat  Plate  Grid  Generation 

GBL  incorporates  two  sample  cases  to  test  the  coding  accuracy  and/or  the  general 
behavior  of  present  and  future  algorithms.  The  first  test  case  simulates  self-similar 
Falkner-Skan  flow  while  the  second  case  simulates  a  linearly  decelerated  flow  with  mild 
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separation,  initially  computed  by  Klineberg  and  Steger.  Both  test  cases,  which  are 
explained  in  detail  later  in  this  report,  require  a  flat  plate  grid  generated  internally  by 
GBL  according  to  user  specifications. 

The  run  parameter  file  contains  the  following  flat  plate  grid  information:  start  coordinate 
ss>  end  coordinate  se,  the  total  number  of  streamwise  stations  NXS,  and  the  number  of 
points  normal  to  the  streamwise  direction  at  each  station  KBL.  GBL  assumes  that  the 
stagnation  point  is  located  at  the  ?  =  0  position  and  that  ss  and  se  are  positive  numbers. 

At  each  station,  the  first  point  above  the  flat  plate  surface,  h2,is  located  at  n+  =  1,  based 
on  the  Prandlt  logarithmic  law  of  the  wall,  and  multiplied  by  the  user  selected  factor  //. 
Freestream  conditions  and  the  theoretical  Blasius  wall  shear  stress  distribution  are  used 
in  the  computation.  Different  values  are  obtained  for  laminar  and  turbulent  flow  as 
follows: 


n2  = 


R  e-frj* 


(2-8a) 


and  the  non  dimensional  wall  shear  stress  Xw  takes  different  values  for  laminar  and 
turbulent  flow.  For  laminar  flow,  it  is  defined  as 


xw/  =  0.332  Ml 

while  for  turbulent  flow,  it  is  given  by 


ujsl 


1/2 


(2-8b) 


xW(  =  0.0225  Mi 


85  *uc 


1/4 


(2-8c) 


where  5*  is  the  turbulent  flat  plate  boundary  layer  displacement  thickness  (see  reference 
10,  p.599) 


5  =0.125  x 


*•  ^  ■> 
Moo?/ 

-1/5* 

0.37?/ 

Voo 

k  J 

(2-9) 


For  each  station,  the  upper  boundary  in  the  n  direction,  is  located  at  the 

corresponding  Blasius  displacement  thickness  multiplied  by  an  amplification  factor  fu  of 
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the  order  of  1.5  to  5.  For  laminar  flow,  this  yields 

*KBL  =fu  5.0 


-O  1/2 
v„  si 


Me 


while  for  turbulent  flow,  we  have 


hkbl  =  L  8* 


(2- 10a) 


(2- 10b) 


The  location  of  the  KBL-3  points  between  «2  an(*  kkbl  ls  determined  from  the 
exponential  stretching  function 

**  =  «*_!  +  h2(\  +£)(k~2)  k  =  3,KBL-l  (2-11) 

where  e  is  computed  using  a  Newton-Raphson  method  to  yield  nKBL  when  k  =  KBL. 
Figure  1  shows  a  typical  flat  plate  grid  generated  by  GBL. 


Figure  1.  Flat  plate  grid  generated  by  GBL 


UNCLASSIFIED 


10 


UNCLASSIFIED 


Navier-Stokes  and  Boundary  Layer  Grids 

When  computing  airfoil  flow,  GBL  requires  two  files  from  the  ARC2D  Navier-Stokes 
code:  the  grid  file  which  contains  the  Cartesian  coordinates  (x,y)  of  the  numerical  grid, 
and  the  state  vector  file  which  contains  the  density,  Jc-momentum,  y-momentum,  and 
energy  information  for  each  grid  point. 

Only  grids  of  C-type  topology,  as  shown  in  figure  2,  are  used  in  the  present  study.  These 
are  generated  with  a  hyperbolic  grid  generation  program.  Each  r\  line  starts  from  the 
lower  far-field  wake,  wraps  around  the  airfoil,  and  ends  at  the  upper  far  wake  to  form  a 
"C"  shape.  The  %  lines,  perpendicular  to  the  rj  lines,  extend  from  the  wake  centerline,  or 
airfoil  surface,  to  the  far  field  outer  boundary.  The  line  defining  the  airfoil  surface  and 
the  wake  centerline  is  the  T]  =  1  line.  The  1  line  is  the  lower  far  wake  boundary.  An 
effort  is  made  to  keep  the  orthogonality  of  each  (£,r|)  line  intersection.  This  is 
particularly  important  near  a  solid  boundary  or  at  the  wake  centerline  where  we  want  to 
solve  the  boundary  layer  equations.  The  reader  should  also  notice  that,  for  a  Navier- 
Stokes  grid,  the  tj  grid  lines  are  concentrated  where  viscous  effects  dominate  the  flow, 
i.e.  near  the  airfoil  surface  and  wake  centerline. 
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Figure  2.  Navier-Stokes  grid  in  Cartesian  coordinates 


Compared  to  the  chord  length,  the  thickness  of  an  airfoil  boundary  layer  is  thin,  growing 
from  a  fraction  of  a  percent  near  the  leading  edge  stagnation  point,  to  some  five  or  ten 
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percent  at  the  trailing  edge.  It  is  then  reasonable  to  assume  that  points  near  the  airfoil 
surface,  for  a  given  ^  station,  are  lined  up  on  the  normal  to  the  surface  at  that  same 
station.  A  close  examination  of  the  Navier-Stokes  grid  shows  this  assumption  to  be  valid 
over  most  of  the  airfoil  except  for  the  trailing  edge  (see  figure  3).  At  that  point,  three  line 
segments  intersect:  two  from  the  airfoil  lower  and  upper  surfaces  respectively,  and  one 
from  the  wake  centerline.  This  causes  a  significant  slope  discontinuity  which  makes  it 
impossible  to  define  a  normal  to  both  the  airfoil  surface  and  the  wake  centerline;  an 
average  normal  is  used  instead.  Hence,  one  must  be  careful  in  using  trailing  edge 
Navier-Stokes  coordinates  for  the  boundary  layer  solution,  since  they  are  not  truly 
normal  to  the  surface.  Furthermore,  the  £,  lines  are  normally  clustered  at  the  trailing  edge 
to  resolve  local  velocity  gradients,  thus  several  lines  may  not  be  normal  to  the  surface  in 
that  region. 


GBL  uses  the  ARC2D  Cartesian  coordinates  ( x,y )  to  generate  one  of  two  possible 
boundary  layer  grids.  The  first  grid  type  consists  of  a  subset  from  the  Navier-Stokes  grid 
while  the  second  grid  type,  discussed  later  in  this  section,  is  an  "adaptive"  grid. 


Generating  the  boundary  layer  grid  (see  figure  4)  from  the  Navier-Stokes  grid  presents 
the  advantage  that  the  Navier-Stokes  and  boundary  layer  procedures  share  the  same  grid. 
Interpolation  of  the  variables  is  avoided  and  the  overall  computing  expense  of  the  FNS 
procedure  is  reduced.  However,  to  be  useful  as  a  boundary  layer  grid,  a  subset  of  the 
Navier-Stokes  grid  grid  must  have  a  sufficient  concentration  of  grid  lines  near  the  airfoil 
surface  to  resolve  the  strong  normal  gradients  characterizing  the  boundary  layer  region. 
This  is  particularly  important  near  the  leading  edge  where  the  boundary  layer  is  very 
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thin.  On  the  other  hand,  typical  Navier-Stokes  grids  are  generated  such  that  the  normal 
distance  between  each  rj  line  and  the  airfoil  surface  is  approximately  constant.  This 
means  only  a  small  number  of  points  span  the  boundary  layer  near  the  leading  edge, 
while  more  points  than  required  span  the  boundary  layer  at  the  trailing  edge.  At  least  ten 
points  should  span  the  boundary  layer  near  the  leading  edge  while  no  more  than  35  to  45 
points  are  required  for  the  trailing  edge.  Fortunately,  the  rapid  stretching  of  the  normal 
coordinates  used  to  generate  the  Navier-Stokes  grid  respects  the  above  conditions  on  the 
number  of  points  across  the  boundary  layer. 

The  boundary  layer  grid  of  figure  4  is  used  only  partially  by  GBL.  For  each  E,  station, 
GBL  determines  the  approximate  number  of  grid  points  that  span  the  boundary  layer  and 
stores  it  in  the  grid  index  array  NDT().  This  is  done  for  two  reasons.  First,  why  waste 
computer  power  on  grid  points  which  are  clearly  outside  the  boundary  layer?  There  is 
however  a  glitch  with  this  approach.  The  finite  difference  algorithm  requires  the 
streamwise  velocity  from  adjacent  stations  to  compute  the  convection  derivatives,  and  at 
times,  these  values  are  from  points  outside  the  boundary  layer.  It  is  then  necessary  to 
update  the  value  of  points  outside  the  NDT()  boundary.  For  the  momentum  and  energy 
equations,  outside  values  are  set  equal  to  the  value  at  NDT().  In  most  cases,  this  is 
acceptable  because  the  NDT()  boundary  varies  smoothly  from  station  to  station. 


Figure  4.  Boundary  layer  grid,  in  Cartesian  coordinates,  formed  from  a  subset  of  the 
Navier-Stokes  grid. 


The  second  reason  to  restrain  the  computations  to  NDT()  points  only  is  related  to  the 
continuity  equation.  The  boundary  layer  equations  are  not  applicable  to  the  flow  outside 
the  boundary  layers  and  wake,  where  they  predict  constant  streamwise  velocity  gradients 
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in  the  outer  flow.  Integration  of  these  gradients  with  the  continuity  equation  yields 
unbounded  normal  velocity  t>  and  leads  to  instability  of  the  numerical  algorithm.  It  is 
then  necessary  to  neglect  the  values  of  v  outside  the  NDT()  boundary. 


Adaptive  Boundary  Layer  Grid 

The  previous  grid  has  only  a  small  number  of  grid  points  in  the  leading  edge  region. 
This  may  affect  the  accuracy  in  the  upstream  region  and  the  resulting  errors  are 
propagated  downstream  due  to  the  parabolic  nature  of  the  equations.  An  "adaptive"  grid, 
where  the  word  adaptive  indicates  that  scaling  of  the  grid  is  such  that  an  approximately 
constant  number  of  points  span  the  physical  boundary  layer  at  any  station,  may  increase 
the  accuracy  in  the  leading  edge  region  and  improve  the  overall  accuracy  of  the  code. 

Several  schemes  were  studied  to  generate  an  adaptive  boundary  layer  grid.  The  one 
selected  uses  the  Navier-Stokes  grid  as  an  mderlay.  This  ensures  the  orthogonality 
characteristics  of  the  Navier-Stokes  grid  are  retained  in  the  boundary  layer  grid,  an 
important  fact  near  the  trailing  edge.  The  overlay  property  also  reduces  the  interpolation 
of  the  Navier-Stokes  data  onto  the  boundary  layer  points  to  a  single  dimension.  An 
example  grid  is  shown  in  figure  5. 
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Figure  5.  Adaptive  boundary  layer  grid  in  Cartesian  coordinates 

The  upper  limit  of  the  adaptive  grid,  in  the  direction  normal  to  the  surface,  is  generated 
from  the  turbulent  flat  plate  boundary  layer  thickness  distribution  of  equation  (2-9)  and 
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Demonstration  grid  for  NACA  0012 
airfoil  at  1.49°  angle  of  attack.  The 
stagnation  point  is  slightly  below  the 
leading  edge. 
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Figure  6.  Adaptive  grid  details  near  the  leading  edge 


Demonstration  grid  for  NACA  0012 
airfoil  at  1.49°  angle  of  attack.  The 
grid  overlays  that  from  ARC2D. 
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Figure  7.  Adaptive  grid  details  near  the  trailing  edge 
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the  appropriate  freestream  conditions.  This  is  adequate  to  approximate  the  growth  of  the 
boundary  layer  over  the  airfoil.  Recall  that  a  flat  plate  experiences  zero  pressure  gradient 
throughout  while  an  airfoil  experiences  a  strong  negative  pressure  gradient  over  its 
forward  portion  followed  by  a  positive  pressure  gradient,  sually  starting  near  be  point 
of  maximum  thickness.  Therefore,  the  flat  plate  distribution  overestimates  the  growth  of 
the  boundary  layer  near  the  airfoil  leading  edge  while  the  opposite  is  true  aft  of  the  point 
of  maximum  thickness.  Hence,  to  ensure  that  the  grid  spans  the  full  boundary  layer,  it  is 
necessary  to  multiply  the  flat  plate  distribution  by  the  factor  fu  which  assumes  a  value 
between  1.5  and  5  depending  on  the  angle  of  attack.  It  should  also  be  noted  that  the 
variable  s  is  defined  as  the  arc  length,  measured  from  the  upstream  stagnation  point, 
along  the  airfoil  surface  or  wake  centerline.  In  the  wake  region,  the  value  of  s  gets  larger 
than  is  intended  with  the  formula  of  equation  (2-9).  Use  of  the  square  root  of  s  results  in 
a  more  realistic  approximation  to  the  growth  of  the  wake. 

The  first  grid  point  above  the  surface  is  set  at  a  constant  distance  «niin  in  the  range 
lxlO-5  to  lxlO-7.  This  number  is  based  on  the  assumption  of  unit  chord.  From 
experience  gained  with  Navier-Stokes  codes,  it  is  known  that  locating  the  first  point  in 
this  range  results  in  computed  value  of  n+  near  or  less  than  one.  In  the  present  boundary 
layer  code,  the  user  specifies  the  distance  of  the  first  grid  point  above  the  airfoil  surface. 
This  value  remains  constant  for  all  streamwise  stations.  Since  the  upper  limit  at  each 
station  grows  as  one  moves  downstream,  while  the  number  of  points  across  the  boundary 
layer  remains  constant,  the  stretching  factor  of  equation  (2-11)  to  locate  the  KBL-3 
points  between  n  2  and  is  recomputed  at  each  station.  The  resulting  values  are  used 
to  interpolate  the  (x ,y)  coordinates  of  the  boundary  layer  grid  from  the  Navier-Stokes 
grid. 

Special  treatment  of  the  stagnation  point  is  required  because  equation  (2-9)  is  singular 
when  j  is  zero.  The  upper  grid  limit  is  then  computed  from  a  linear  extrapolation  of  the 
neighbouring  points.  Figure  6  shows  the  result.  The  grid  lines  obtained  near  the  trailing 
edge  are  shown  in  figure  7. 

Grid  generation  schemes  other  than  the  one  above  were  studied,  but  were  not  found  to  be 
as  successful.  For  example,  computing  the  upper  boundary  distance  by  using  the 
Navier-Stokes  data  to  estimate  local  boundary  layer  thickness  does  not  work  well 
because  there  is  no  clear  distinction  between  the  boundary  layer  and  the  inviscid  flow 
region.  An  attempt  to  compute  the  location  of  the  first  grid  point  above  the  surface  using 
the  wall  shear  stress  was  also  made,  but  was  not  found  to  provide  reliable  results  near  the 
stagnation  point  and  in  regions  of  separated  flow.  It  was  therefore  concluded  that, 
although  not  optimal,  the  above  grid  generation  scheme  yields  a  smooth  grid  and  does  it 
efficiently. 

After  generating  the  Cartesian  grid  coordinates  ( x,y ),  GBL  transforms  them  to  body 
conformal  coordinates  (s,n)  as  required  by  the  current  boundary  layer  scheme.  The  s 
coordinates  are  obtained  from  integration  of  the  T|-  1  lines.  For  any  given  £,  station,  the 
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value  of  s  remains  constant.  The  s  =  0  coordinate  is  located  at  the  leading  edge 
stagnation  point  with  the  positive  coordinates  corresponding  to  the  upper  surface  of  the 
airfoil.  The  n  coordinates,  defined  to  be  zero  on  the  airfoil  surface  or  wake  centerline, 
are  obtained  by  integrating  each  4  line  from  the  TJ  =  0  point. 

The  integration  process  to  generate  the  body  conformal  coordinates  is  equivalent  to 
unwrapping  the  Cartesian  grid  such  that  the  T]  =  1  line  is  straight.  It  also  means  that  T| 
lines  above  the  surface  are  compressed  (their  total  arc  length  is  reduced)  in  the  process, 
but  this  is  neglected  because  the  boundary  layer  normal  dimensions  are  very  small  when 
compared  to  the  streamwise  dimensions.  A  typical  body  conformal  grid  is  shown  in 
figure  8. 


s  coordinate 


Figure  8.  Boundary  Layer  grid  in  body  conformal  coordinates 
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COMPUTATIONAL  TOOLS 

Several  computational  tools  are  essential  to  the  success  of  the  boundary  layer  algorithm. 
For  example,  one  tool  is  required  to  compute  the  metrics  of  the  transformation  when 
mapping  the  physical  grid  to  the  computational  grid.  Other  tools  transform  the  velocity 
components  from  one  coordinate  system  to  another,  compute  vorticity,  or  implement  the 
Baldwin-Lomax  turbulence  model^1^  (suitably  modified  for  the  body  conformal 
coordinate  system  used  in  GBL).  The  main  algorithm  requires  the  solution  of  linear 
systems  of  equations  at  each  streamwise  station;  a  tridiagonal  solver  is  used  for  the  direct 
mode  while  the  inverse  mode  requires  a  modified  version  of  the  same  algorithm.  Finally, 
other  tools  yield  an  estimate  of  the  convergence  of  the  algorithm. 


Metrics  of  the  Transformation 

Use  of  the  Cartesian  (x,y)  or  body  conformal  (s,h)  grid  coordinates  to  solve  the  boundary 
layer  equations  increases  the  complexity  of  the  finite  difference  algorithm  because  the 
spacing  between  points  is  not  uniform.  It  is  preferable  to  transform  the  equations  to 
solve  them  on  a  uniform  grid  with  unit  spacing  in  both  directions.  This  increases  the 
complexity  of  the  equations  through  the  addition  of  metric  terms,  but  simplifies  the  finite 
difference  approximations  to  the  derivatives  and  results  in  an  overall  gain  in 
computational  efficiency. 

The  physical  grid  points  are  mapped  one-to-one  to  the  computational  grid  with 
coordinates  (£,,ri).  The  only  exception,  for  the  C-grid  used  in  the  present  study,  is  the 
wake  cut  where  each  physical  point  is  mapped  to  two  computational  points.  Conversely, 
one  can  think  of  it  as  two  computational  points  that  overlap  in  the  physical  grid.  Figure  9 
shows  an  example  of  the  physical  (body  conformal  coordinates)  and  computational  grids. 


Body-Conformal 


Computational 


Figure  9.  Body  conformal  and  computational  grids 

As  mentionned  above,  the  transformed  equations  contain  metrics  of  the  transformation 
(^x,  Tlx  and  Tl>)  which  must  be  determined  numerically.  This  would  normally  involve 
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finite  differences  using  the  physical  grid  for  which  the  £  and  T)  coordinates  corresponding 
to  each  point  are  known.  However,  these  computations  use  unequal  grid  spacing  once 
again,  exactly  what  we  are  trying  to  avoid.  It  is  more  practical  to  approach  the  problem 
from  the  other  end.  The  physical  coordinates  corresponding  to  each  (4,q)  line  intercept 
in  the  computational  domain  are  known,  so  why  not  compute  the  derivatives  x^,  x^, 
and  y^  using  the  computational  grid,  and  then  compute  the  metrics  from  these?  All  that 
is  required  is  the  following  relation  linking  the  computational  plane  derivatives  to  the 
metrics. 

The  functional  form  of  the  computational  variables  is 


S  =  k$,y) 

T1  =  Tlfty) 


(3-1) 


Using  the  chain  rule,  differential  increments  in  the  computational  domain,  d£,  and  dq,  are 
linked  to  increments  dx  and  dy  in  the  physical  space,  i.e. 


Si  Sjl  [t& 

"Hi  dy 


(3-2) 


The  role  of  the  dependent  and  independent  variables  may  also  be  reversed,  i.e. 


x=x(^,q) 

5;=5'(tTl) 

for  which  the  differential  increments  relation  is 


'dx\  *n  \d% 

dy  y^  y^  dr> 


(3-3) 


(3-4) 


Solving  this  last  equation  for  d ^  and  dq  in  terms  of  dx  and  dy,  and  equating  with 
equation  (3-2),  results  in  the  relation 
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'fe  4>' 

which,  in  its  expanded  form  gives  the  metrics 

4i=^“15'  r\  .  T\i=-J~xh 

ky=-J~xxr[  ,  r\y=  J~xx^ 

where  J,  the  Jacobian  of  the  transformation,  is  defined  as 


*4 

h  h 


(3-5) 


(3-6) 


•/  =  (jv5 


Cartesian  to  Body  Conformal  Velocity  Transformation 

During  the  initialization  process  for  airfoil  computations,  GBL  reads  the  ARC2D  file 
inputs.  These  are  in  Cartesian  coordinates,  while  GBL  uses  body  conformal  coordinates 
for  its  computations.  Scalar  quantities,  such  as  state  and  energy  variables,  are 
independent  of  the  coordinate  system  used,  but  vector  quantities  (such  as  velocity  or 
momentum)  are  not  and  must  be  transformed  from  one  coordinate  system  to  another. 
The  transformation  requires  the  cosine  angles  relating  the  two  coordinate  systems.  They 
are  defined  from  the  assumption  that  during  grid  generation,  points  at  a  given  E,  station 
were  located  on  the  normal  to  the  airfoil  surface  (see  Figure  10),  with  the  minor 
exception  of  a  few  lines  in  the  trailing  edge  region  as  discussed  previously. 
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The  metrics  of  the  transformation  rjj  and  T]y,  mapping  the  Cartesian  grid  coordinates  to 
the  computational  domain,  are  then  the  x  and  y  components  of  the  vector  q  normal  to  the 
airfoil  surface.  The  angle  0  represents  the  rotation  angle  between  the  Cartesian  and  body 
conformal  coordinate  systems.  Using  the  metrics,  the  sine  and  cosine  of  0  are  defined  as 


COS0  = 


sin0  = 


15 


(Tljl+Tl?)* 

-m 

Oli+Tly)* 


(3-7) 


These  trigonometric  quantities  relate  the  Cartesian  and  body  conformal  components  of 
the  velocity  vector  through  the  equations 


ub  =  uc  cos0  +  sin0 
\j b=-uc  sin0  +  cos0 


(3-8) 


The  inverse  relations  are 


uc  =  ub  cos0  -  vb  sin0 
oc  =  ub  sin0  +  vb  cos0 


(3-9) 


The  transformation  sine  and  cosine  angles  are  computed  during  initialization  and  stored 
in  arrays  SINANGO  and  COSANG().  The  transformation  from  Cartesian  to  body 
conformal  coordinates  is  handled  by  subroutine  UVTC2B  while  the  inverse  transform, 
used  during  output,  is  done  by  subroutine  UVTB2C. 


Vorticity  Computations 

During  initialization  of  airfoil  cases,  the  vorticity  is  used  to  interpolate  the  ARC2D  data 
in  the  q  direction,  for  each  %  station,  to  find  the  location  where  the  magnitude  of  the  local 
vorticity  w  is  less  than  or  equal  to  the  user  supplied  cutoff  value  I  tbcu/ 1 .  The  vorticity  is 
defined  in  terms  of  Cartesian  velocity  components.  In  its  non-dimensionalized  form,  it  is 


dv 

to  =  — 

dx 


du 

dy 


(3-10) 
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The  metrics  of  the  transformation  from  Cartesian  to  computational  coordinates  are  used 
to  replace  the  x  and  y  derivatives  and  result  in  the  expression 

©  =  0.5^5+11^  -^yU^-T]yU^  (3-11) 

which  is  evaluated  with  centered  differences  wherever  possible.  At  the  boundaries, 
backward  or  forward  three  point  differences  are  used. 


Turbulence  Modelling 

The  turbulent  flow  equations  are  derived  from  the  laminar  flow  equations  by  replacing 
velocity,  pressure,  density  and  temperature  by  the  sum  of  their  mean  and  fluctuating 
components,  i.e.  u  =  u+u p=p  +p',  etc.  A  time  average  of  the  equations  is  then 
taken  and  negligible  terms  are  dropped.  The  resulting  equations  are  identical  to  those  for 
laminar  flow  except  for  additional  terms  such  as  the  Reynolds  stresses  -pu'v  .  These 
express  the  increased  diffusion  characteristics  of  turbulent  flow. 

The  computational  algorithm  of  GBL  applies  to  laminar  as  well  as  turbulent  flow 
provided  the  increased  diffusion  is  considered.  For  the  momentum  equation,  the 
viscosity  p  is  replaced  by  the  sum  of  the  molecular  viscosity  pm,  obtained  from 
Sutherland’s  formula,  and  a  turbulent  viscosity  p,  which  replaces  the  Reynolds  stresses. 
In  the  energy  equation,  the  ratio  pJPr  is  replaced  by  )AmIPrm  +  ptIPrt  where  Prm  =  0.72 
and  Pr,  =  0.9  are  the  laminar  and  turbulent  Prandtl  numbers,  respectively.  A  turbulence 
model  is  required  to  evaluate  p,. 

The  Baldwin  and  Lomax* 1 !  *  model  is  widely  used  with  the  thin-layer  Navier-Stokes 
equations.  It  uses  a  two  layer  algebraic  eddy  viscosity  model  to  evaluate  p,.  In  the  inner 
region,  near  a  solid  wall,  it  uses  the  Prandtl-Van  Driest  formulation.  In  the  outer  region, 
and  in  wakes,  it  uses  the  Clauser  formulation.  Transition  from  one  formulation  to  the 
other  is  done  at  a  crossover  point,  above  the  wall,  defined  as  the  first  point  where  the  two 
formulations  are  equal.  The  reader  is  referred  to  the  original  paper  for  the  development 
and  justification  of  the  method. 

The  original  Baldwin-Lomax  equations  are  expressed  below  in  terms  of  the  scaled 
variables  used  in  GBL.  The  inner  region  formulation  is  given  by 

p{(.  =  Re«- p -d2-  ld>l  (3- 1 2a) 

The  variable  d  is  a  characteristic  length  scale  of  the  turbulence  defined  as 
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cf  =  0.4n  (1  -  e-n+/26) 


(3-l2b) 


where  n+  is  a  local  Reynolds  number  defined  as 


n+  =h 


Pw^w 

Re  <,,—72 


>/4 


(3- 12c) 


In  the  outer  region,  the  scaled  turbulent  viscosity  is  given  by 


PtQ  -  0.0168-  1.6-  p  -  Fv/ake'pKleb 


(3- 13a) 


The  function  F/oebin)  is  the  Klebanoff  intermittency  factor  defined  as 


1+5.5 

0.3  n 

6" 

_ 

.  nmax  _ 

_ 

(3-1 3b) 


while  the  function  Fwak£  is  defined  as 


Fwake  =  MIN\ 


ft  max  ^  max  > 


ftma x^dif 


(3-13c) 


The  quantities  nmax  and  Fmax  are  determined  from  the  function 

F(n)  =  n  ! co I  ( 1  -  e~n+l 26) 


(3- 13d) 


and  the  quantity  udlf  is  the  difference  between  the  maximum  and  minimum  total  velocity 
in  the  profile,  i.e. 


„2  _2  Vi  _  2  _2  V4 

Udif  ~  (W  +  V)  Jfnajt  —  (u  +  \)  )min 


(3-130 
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Interpolation,  using  a  local  second  order  polynomial  fit,  improves  the  values  of  nmax  and 
Fmax.  1°  particular,  the  variation  of  Fmax  is  smoother  in  going  from  one  station  to  the 
next. 

The  above  model  assumes  Cartesian  coordinates.  It  is  modified  for  use  with  body 
conformal  coordinates.  First,  the  uc  and  \)c  components  of  velocity  are  replaced  by  their 
body  conformal  counterparts  ub  and  X)b-  Furthermore,  outside  the  boundary  layer  edge 
defined  by  the  index  NDT,  \>b  is  neglected  because  it  is  unbounded  outside  this  level,  as 
explained  before.  This  is  particularly  true  near  the  leading  edge  region  where  the  ub 
velocity  gradient  is  very  strong  and  the  boundary  layer  is  thin. 

The  second  change  is  the  replacement  of  the  vorticity  by  the  normal  gradient  of  ub,  i.e. 
I  to  I  — )  \du  /  dn  I  The  behavior  of  these  two  quantities  is  similar  with  differences  in  the 
magnitude  of  the  terms,  but  use  of  the  body  conformal  components  compensates  for  that. 


Direct  and  Inverse  Solvers 

The  direct  mode  of  GBL  uses  a  finite  difference  scheme  with  centered  or  three  point 
forward/backward  difference  stencils.  The  equations  are  solved  at  each  t,  station,  which 
results  in  a  scalar  tridiagonal  systems  of  equations  in  the  rj  direction,  i.e. 

lA$=t 

where  [A]  is  a  tridiagonal  matrix,  it  is  the  vector  of  unknowns  to  be  solved  for,  and  li  is 
the  right  hand  side  vector,  a  known  quantity.  The  efficient  Thomas  algorithm112*  is  used 
to  establish  an  upper  triangular  matrix  followed  by  back  substitution  to  yield  a  solution. 
Subroutine  THOMAS,  presented  in  Appendix  B,  solves  the  above  system  of  equations. 

For  thejnverse  mode,  the  pressure  forcing  function,  which  is  buried  in  the  right  hand  side 
vector  B,  must  be  replaced  by  an  inverse  function.  GBL  implements  the  inverse  forcing 
functions  of  Van  Dalsem  and  Steger*5*.  They  apply  the  momentum  equation  at  the  wall  or 
on  the  wake  centerline  to  derive  an  expression  to  replace  the  pressure.  The  resulting 
expression  is  a  function  of  the  wall  shear  stress  over  the  airfoil,  and  of  the  wake 
centerline  velocity  in  the  wake.  This  yields  a  modified  tridiagonal  system  of  equations 
which  has  non -zero  entries  in  the  first  column  as  shown  below.  A  lower-upper 
decomposition  scheme  is  used  to  solve  this  particular  matrix  equation;  back  substitution 
is  applied  to  the  upper  triangular  matrix  to  obtain  an  intermediate  solution  which  is  then 
used  with  the  lower  triangular  matrix  and  forward  substitution  to  yield  the  final  result. 
The  subroutine  listing  is  included  under  the  name  INVTOM  in  Appendix  B. 
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Convergence  Estimates 

It  is  useful  to  have  an  estimate  of  the  convergence  of  the  numerical  algorithm.  GBL 
implements  two  methods  to  produce  such  an  estimate.  The  first  method  takes  a  simple 
difference  of  representative  field  variables  from  one  cycle  to  the  next,  i.e.  variable  at  time 
n  minus  variable  at  time  n— 1.  The  quantity  used  for  each  equation  is:  u  for  momentum, 
H  for  energy,  and  p  for  continuity.  The  process  involves  all  points  of  the  field.  The 
location  and  magnitude  of  the  largest  change  is  saved  and  the  root  mean  square  for  the 
entire  field  is  computed.  Its  logarithm  (base  ten)  is  presented,  i.e. 


A  Var  =  logio(RMS) 


(3-14) 


The  second  method  substitutes  the  variables  back  into  their  finite  difference  equations 
after  each  cycle.  As  the  algorithm  converges,  the  difference  between  the  left  and  right 
hand  sides  of  the  equations  approaches  machine  zero.  Only  interior  points  are  used  and 
equations  (3-14)  and  (3-15)  remain  valid.  Finally,  when  interpreting  results  obtained 
from  either  method,  the  user  must  be  aware  that  convergence  depends  on  the  magnitude 
of  Ar,  the  time  step  used  to  relax  the  equations. 
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INITIALIZATION  OF  VARIABLES 

The  numerical  method  used  in  GBL  requires  initialization  of _ the  forcing  function 
variables  ipe^w  and  uwc),  boundary  condition  variables  ( ue  and  He)  and  field  variables 
( u ,  t),  U ,  V,  H,  p  and  p).  The  algorithm  uses  these  initial  values  to  estimate  flow  field 
derivatives  and  start  the  iterative  cycle.  Furthermore,  because  of  to  the  non-linear  nature 
of  the  boundary  layer  equations,  the  initial  estimates  must  be  representative  of  the  flow 
field  being  computed;  otherwise,  the  method  may  not  converge.  This  type  of  behavior  is 
characteristic  of  non-linear  systems. 

GBL  supports  four  initialization  procedures,  two  for  test  purposes  and  two  for  airfoil 
computations.  The  test  cases  are  for  laminar  flow  only.  They  are  Falkner-Skan  and 
Klineberg  flows.  The  Falkner-Skan  case  is  used  to  check  the  direct  mode  portion  of 
GBL.  The  program  uses  theoretical  distributions  of  pressure  and  velocity  to  initialize  the 
flow  field.  The  Klineberg  case  is  used  to  exercise  the  inverse  mode  portion  of  GBL.  It  is  a 
linearly  retarded,  mildly  separated  flow  for  which  the  distribution  of  wall  shear  stress  is 
provided.  The  flow  field  is  initialized  with  scaled  Blasius  profiles.  It  should  be  noted 
that  although  both  cases  are  valid  for  laminar  flow  only;  the  Falkner-Skan  test  case  can 
also  be  used  to  check  the  implementation  of  the  turbulence  model.  All  test  case  use  the 
flat  plate  grid  generated  internally  by  GBL  (as  described  in  section  2). 

The  initialization  process  for  airfoil  flow  is  compatible  with  the  ARC2D  Navier-Stokes 
code.  The  procedure  uses  the  ARC2D  state  vector  and  the  corresponding  grid  to  generate 
an  "adaptive"  grid  or  uses  a  subset  of  the  Navier-Stokes  grid.  In  the  former  case,  the 
ARC2D  data  is  interpolated  linearly  to  initialize  field  values  at  the  boundary  layer  grid 
points.  For  the  latter  case,  the  ARC2D  field  data  is  used  directly.  A  coordinate 
transformation  is  then  applied  to  the  velocity  data. 


Velocity-Pressure  Relation  at  Boundary  Layer  Edge 

The  Falkner-Skan  test  case  has  ue(s)  prescribed,  but  the  corresponding  edge  pressure 
distribution,  pe(s),  is  unknown.  It  is  the  opposite  for  the  airfoil  initialization  cases.  The 
surface  pressure  from  the  Navier-Stokes  solution  may  be  used  as  pe(s )  if  it  is  assumed 
constant  across  the  boundary  layer.  However,  finding  ue(s )  is  difficult  because  the 
Navier-Stokes  solution  does  not  show  a  clear  transition  between  the  boundary  layer  and 
the  inviscid  flow  field.  A  relationship  between  we(s)  and  pe(s )  is  then  required. 

The  boundary  layer  effectively  displaces  the  streamlines  of  the  inviscid  flow  away  from 
the  solid  surface  by  a  distance  equal  to  the  displacement  thickness  of  the  boundary  layer. 
If  this  fact  is  combined  with  the  observation  that  the  pressure  remains  constant  across  the 
boundary  layer,  Euler’s  inviscid  momentum  equation,  which  is  valid  along  a  streamline, 
may  then  be  used  to  relate  the  edge  velocity  distribution  ue(s )  to  the  surface  pressure 
distribution  pe(s). 
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(4-1) 


The  above  constant  is  computed  by  replacing  the  edge  values  of  velocity,  pressure  and 
density  by  the  freestream  equivalents.  It  assumes  the  flow  is  isentropic,  i.e.  it  satisfies  the 
relation 


p  =  (constant)p (4-2) 

which  holds  true  for  most  subsonic  and  transonic  flows  under  consideration  in  this  study. 
Equation  (4-2)  is  used  to  eliminate  pe  in  favor  of  pe  in  equation  (4-1).  Scaling  the 
variables  in  the  resulting  equation,  we  obtain  expressions  linking  pressure  and  velocity  at 
the  boundary  layer  edge. 


pe  =  (pa£)  1 


p~ 

uL  -  (a^Ue)2 

*  •>. 

y- 1 
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(4-4) 


An  expressions  to  compute  the  temperature  at  the  boundary  layer  edge  is  obtained  by 
combining  the  isentropic  relationship  with  the  perfect  gas  law 


(4-5) 


The  temperature  is  then  used  to  compute  total  enthalpy  at  the  boundary  layer  edge 


Pe  P 


2  ^ 


(4-6) 
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Falkner-Skan  Test  Case 

A  flow  field  is  called  self-similar  when  velocity  profiles  u(s,h)  at  different  stations  s 
differ  only  by  scale  factors  in  u  and  h.  Schlichting*10^  (pp.  152- 156)  states  the  conditions 
under  which  self-similar  conditions  exist.  One  class  of  self-similar  flow,  studied 
extensively  by  Falkner  and  Skan^13\  is  obtained  by  specifying  the  boundary  layer  edge 
velocity  as 


ue(s)  =  Csm ,  m  =  (4-7) 

This  function  describes  the  potential  flow  in  the  neighbourhood  of  the  stagnation  point  of 
a  wedge  whose  included  angle  is  7$.  The  distance  s  along  the  surface,  is  defined  to  be 
zero  at  the  stagnation  point.  A  non-dimensional  form  of  this  equation,  with  the  constant 
C  set  equal  to  «<*,,  is  used  in  GBL. 


ue=M00(sl)m 


(4-8) 


With  the  condition  of  equation  (4-7),  the  continuity  and  streamwise  momentum  equations 
reduce  to  a  single  ordinary  differential  equation  in  terms  of  the  stream  function  /  and  its 
derivatives  with  respect  to  a  non-dimensional  normal  distance  £. 


/'''  +  a//"  +  |3(l-/'2)  =  0 


(4-9) 


subject  to  the  boundary  conditions 


C  =  0  :/  =  /'  =  0  ;  £  =  ° °:/'=l 


(4-10) 


where  f'  =  u/ue 

By  carefully  selecting  values  of  the  constants  a  and  p,  equation  (4-9)  aporoximates 
important  existing  flows.  The  constant  a  is  often  set  equal  to  unity;  this  is  the  only  case 
supported  by  GBL.  If  the  constant  p  is  set  to  zero,  the  equation  describes  flow  over  a  flat 
plate  at  zero  incidence.  The  case  P  =  1  describes  the  flow  in  the  viscinity  of  a  two 
dimensional  stagnation  point.  Selecting  any  other  values  of  P  between  zero  and  one,  the 
equation  describes  the  accelerating  flow  in  the  neighbourhood  of  the  stagnation  point  of  a 
wedge.  Finally,  by  selecting  small  negative  values  of  P,  the  equation  represents 
decelarating  flow;  separation  is  encountered  at  P  = -0.199.  For  values  of  P< -0.199, 
solutions  of  equation  (4-9)  are  no  longer  unique;  they  are  not  attempted  in  GBL. 
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The  user  specifies  the  self-similar  parameter  m  in  the  input  file  (see  section  2).  GBL  uses 
this  value  to  integrate  equation  (4-9)  from  £  =  0  to  £  =  10  (a  useful  approximation  to 
£  =  °o)  with  a  Runge-Kutta  scheme.  The  integration  is  sensitive  to  the  second  derivative 
at  the  wall,  / ",  and  requires  a  shooting  technique  to  update  this  value  until  convergence 
is  achieved.  The  integration  scheme  yields  distributions  for  the  stream  function  /  (£)  and 
its  derivatives  /'( Q  and 

Subroutine  INITO  initializes  the  flow  field  for  Falkner-Skan  test  cases.  The  procedure 
starts  by  generating  the  Cartesian  coordinates  of  the  boundary  layer  grid,  which,  for  this 
case,  are  also  the  body  conformal  grid  coordinates.  The  metrics  of  the  transformation  are 
then  computed  and  equation  (4-9)  is  solved  to  obtain  the  Falkner-Skan  solution  requested 
by  the  user. 

The  boundary  conditions  and  forcing  functions  are  then  computed.  The  edge  velocity 
distribution  ue(s )  is  obtained  from  equation  (4-8)  while  equation  (4-3)  is  used  to  obtain 
the  pressure  pe(s).  Total  enthalpy  is  computed  from  equation  (4-6).  The  wall  shear 
stress  distribution  xw(5)  is  computed  from  the  theoretical  Falkner-Skan  distribution,  i.e. 


m  4-1 

2  Vpo  1 5  I  / 


Vi 


(4-11) 


Finally,  the  wake  centerline  velocity,  which  is  not  used  for  this  case,  is  set  equal  to  zero. 
Having  initialized  the  forcing  function  and  boundary  condition  arrays,  the  velocity 
components  u  and  x>  are  computed  for  each  grid  point.  Using  the  s  and  n  values,  the 
scaled  value  C,  is  computed  from 


c=« 
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(4-12) 


which  is  used  to  interpolate  values  of  /  and  /'.  These  values  are  then  used  in  the 
following  expressions  to  compute  u  and  o 


(4-13) 

(4-14) 
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The  final  steps  of  INITO  are  to  initialize  the  contravariant  velocity  components,  density 
and  viscosity  throughout  the  field.  For  turbulent  flow,  the  turbulent  viscosity  is  also 
computed. 


Klineberg  Test  Case 

Klineberg  and  Steger^14'  used  an  inverse  method  to  solve  the  laminar  boundary  layer 
equations  for  linearly  decelerated  flows  over  a  flat  plate.  The  wall  shear  stress 
distribution  was  specified  as  an  analytic  function  to  generate  a  small  region  of  reversed 
flow.  Their  paper  gives  an  example  wall  shear  stress  distributk 


x 


0.33238 

12 


(J-2X5-6) 


0  <  s  <  2  ,  s  >  6 


x  =  ToO  +0.1(5  -2)(s-6)  ,  2  <5  <6 


(4-15) 


and  the  corresponding  numerical  results. 

Routine  INIT1  initializes  the  flow  field  variables  for  the  Klineberg  case.  The  first 
operation  is  to  generate  the  boundary  layer  grid  and  to  compute  the  metrics  of  the 
transformation.  The  boundary  conditions  and  forcing  functions  are  then  initialized. 
Klineberg’s  edge  velocity  is  read  from  a  file,  rescaled  using  the  freestream  Mach  number, 
and  a  spline  is  fit  through  the  data  to  interpolate  values  corresponding  to  the  boundary 
layer  grid.  Equation  (4-3)  is  then  used  to  compute  the  edge  pressure  and  equation  (4-6) 
for  total  enthalpy  at  the  edge.  The  wall  shear  stress  distribution  of  equation  (4-15)  is 
rescaled  as  follows  for  use  in  GBL 
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(4-16) 


Finally,  the  wake  centerline  velocity,  which  is  not  used,  is  set  equal  to  zero. 

Since  the  Klineberg  case  starts  from  flat  plate  flow,  the  field  velocity  components  u  and  o 
are  initialized  using  flat  plate  profiles  generated  from  the  Faikner-Skan  equation  with 
a  =  1  and  (3  =  0.  The  remaining  variables  are  then  initialized.  The  contravariant  velocity 
arrays  are  computed  from  the  metrics  of  the  transformation  and  the  body  conformal 
velocity  components  u ^  and  X)f,.  Density  is  computed  from  the  perfect  gas  relation, 
viscosity  from  Sutherland’s  law,  and  turbulent  viscosity  is  set  to  zero. 


UNCLASSIFIED 


30 


UNCLASSIFIED 


Airfoil  Cases 

Two  initialization  procedures  were  designed  for  airfoil  cases.  INIT2  uses  a  subset  of  the 
Navier-Stokes  grid  as  the  boundary  layer  grid,  while  INIT3  generates  an  adaptive 
boundary  layer  grid  as  discussed  in  section  2.  The  two  procedures  are  identical  with  the 
following  exception:  the  Navier-Stokes  data  is  used  directly  in  ENIT2,  but  the  use  of  a 
different  boundary  layer  grid  in  INIT3  requires  interpolation  of  the  Navier-Stokes  data. 

The  initialization  process  starts  with  the  input  of  the  Navier-Stokes  grid  Cartesian 
coordinates  and  the  corresponding  state  vector  Q  file  which  contains  the  non-dimensional 
density,  ^-momentum,  y-momentum  and  total  energy  at  each  point,  i.e. 


Q  = 


P 

P  Uc 
~PX>C 
e 


(4-17) 


The  energy  e  is  related  to  total  enthalpy  by  the  relation 


H  =  —  ~  (Y~  1) 
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f  ~2  - 2>) 

1  uc  +  uc  1 


(4-18) 


The  subscript  c  on  the  momentum  elements  of  Q  denotes  that  Cartesian  components  of 
velocity  are  used.  They  are  used  with  the  metrics  of  the  transformation  to  go  from  the 
Navier-Stokes  to  computational  grids,  and  equation  (3-11),  to  compute  vorticity.  The 
metrics  are  then  used  with  equation  (3-7)  to  compute  the  rotation  angles  necessary  to 
transform  the  the  velocity  components  to  body  conformal  coordinates,  and  the  velocity 
components  are  transformed. 

The  magnitude  of  vorticity  is  strongest  near  the  airfoil  surface  and  vanishes  rapidly  as  we 
move  out  of  the  boundary  layer.  It  is  therefore  used  to  estimate  the  location  of  the 
boundary  layer  edge.  The  user  selects  the  cutoff  vorticity  level  1 0dcui  I  which  defines  the 
edge  and  GBL  interpolates  the  vorticity  data  to  find  this  line.  The  corresponding  edge 
pressure  pe  and  velocity  ue  may  also  be  used  as  boundary  conditions  if  the  user  sets  the 
input  variables  RSPEDG  and  RSUEDG  equal  to  1  in  the  input  file.  The  user  may  also 
elect  to  use  the  surface  pressure  and  the  corresponding  edge  velocity  computed  from 
equation  (4-4)  as  edge  conditions.  The  INIT2  routine  may  also  use  the  pressure  and 
velocity  at  ti  =  KBL  as  boundary  conditions;  INIT3  does  not  offer  this  option. 
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The  stagnation  point  is  located  at  the  £  station  with  the  highest  pressure.  In  addition,  the 
Cartesian  grid  data  is  searched  to  identify  the  ^  stations  corresponding  to  the  leading 
edge  and  lower/upper  trailing  edges.  The  origin  of  the  body  conformal  grid  is  then 
moved  to  the  stagnation  point,  so  that  I  s  1  becomes  the  arc  length  distance  from  the 
stagnation  point,  and  the  metrics  of  the  transformation  from  body  conformal  to 
computational  coordinates  are  computed.  The  contravariant  velocity  components  are 
also  computed  at  this  stage. 

The  next  steps  for  the  initialization  are  to  reset  the  total  enthalpy  if  the  constant  enthalpy 
option  is  selected,  compute  the  molecular  and  turbulent  viscosities,  and  if  the  user  elects, 
to  reset  the  velocity  profiles  at  the  stagnation  and  neighbouring  points  to  stagnation 
(m  =  1)  Falkner-Skan  profiles.  This  last  option  is  included  in  GBL  to  study  the  impact  of 
the  upstream  conditions  on  the  computations. 

Since  the  user  selects  the  boundary  layer  edge  arbitrarily  with  tne  variable  l(5cu/l, 
additional  grid  points  are  added  above  the  previously  computed  boundary  layer  edge  to 
make  sure  that  the  grid  spans  the  entire  boundary  layer.  The  wake  centerline  velocity  is 
then  set  equal  to  the  ii(,  velocity  along  the  line  T)  =  1  and  the  edge  total  enthalpy  is  set  to 
the  enthalpy  at  the  points  rj  =  NDT  (£).  The  wall  shear  stress  forcing  function  is  set  to 
zero  in  the  wake  while,  on  the  airfoil  surface,  it  is  computed  using  a  three  point  forward 
difference  of  the  velocity  component  tangential  to  the  wall,  i.e. 


Xw  —  T|y  Aq 


3  -  £ 


+1 


Ub 


(4-19) 
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SOLUTION  SCHEME 

The  momentum,  energy  and  continuity  equations  share  a  common  solution  process.  The 
computational  field  is  swept  in  a  £,  direction,  and  at  each  station,  an  implicit  system  of 
linear  equations  in  the  T)  direction  is  solved.  GBL  takes  advantage  of  this  commonality 
to  establish  three  levels  of  subroutines  for  the  solution  process.  The  highest  level 
controls  the  overall  solution  cycle  by  calling  each  equation  solver  in  a  particular  order. 
The  second  level  solves  a  given  equation  (e.g.  momentum  )  over  portions  of,  or  over  the 
entire  computational  field.  Check  functions  and  forcing  function  updates  are  also  done  at 
the  second  level.  The  third  level  evaluates  specific  functions  for  the  equation  solvers. 
The  above  structure  is  illustrated  in  Figure  11.  It  provides  the  flexibility  required  for  the 
research  framework  of  GBL.  Changes  to  the  formulation,  boundary  conditions  or  overall 
sweeping  schemes  are  easily  implemented  by  modifying  the  appropriate  modules. 
Elements  for  each  level  are  now  discussed. 


Figure  11.  Levels  in  solution  process 


Sweep  Control  Level 

The  solution  scheme  developed  by  Van  Dalsem  and  Steger  involves  the  solution  of  the 
equations  in  a  particular  order  over  a  user  specified  number  of  cycles.  During  each  cycle, 
the  solution  progresses  as  follows: 
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[1]  The  streamwise  velocity  component  u  is  updated  for  a  portion  or  the  full 
computational  field  by  solving  the  momentum  equation.  The  other  equations  are 
solved  for  the  same  extent  of  the  field. 

[2]  The  contravariant  components  of  velocity,  U  and  V,  are  updated  using  their 
definition  and  the  latest  values  of  u. 

[3]  Total  enthalpy,  H,  is  updated  using  the  energy  equation.  Alternatively,  the 
assumption  of  constant  enthalpy  may  be  used. 

[4]  Density,  p,  is  updated  using  the  latest  values  of  total  enthalpy  and  streamwise 
velocity. 

[5]  The  normal  velocity  component  x>  is  updated  by  integrating  the  continuity 
equation  from  the  rj  =  1  line  and  using  the  latest  values  of  the  streamwise  velocity 
components  and  density. 

[6]  The  contravariant  velocity  component  V  is  once  more  updated  using  the  latest 
value  of  v.  U  does  not  depend  on  x>  and  needs  no  updating. 

[7]  Finally,  the  molecular  viscosity  is  computed  from  Sutherland’s  law  while  the 
Baldwin-Lomax  turbulence  model  is  used  to  update  turbulent  viscosity. 


It  is  not  clear  how  the  computational  field  should  be  swept.  The  parabolic  nature  of  the 
boundary  layer  equations  dictates  that  sweeps  be  done  in  a  downstream  direction,  but 
how  are  embedded  regions  of  reverse  flow  handled?  The  current  sweeping  schemes  of 
GBL  attempt  to  answer  this  question.  Three  sweep  procedures  are  incorporated  in  the 
program.  One  routine,  CMPGBO,  is  used  with  the  flat  plate  grid  to  test  algorithms.  The 
user  controls  the  direction  and  range  of  the  sweeps.  The  other  routines  are  used  for 
airfoil  cases. 

Routine  CMPGB1 A  sweeps  a  user  selected  range  of  the  computational  field,  from  one 
end  to  the  other,  in  alternating  directions  from  one  cycle  to  the  next.  The  decision  to 
alternate  sweep  direction  is  based  on  observations  of  the  behaviour  of  the  algorithm  for 
flat  plate  flows.  First,  the  algorithm  converges  when  sweeping  the  flow  field  in  the 
general  downstream  direction,  which  is  consistent  with  the  parabolic  nature  of  the 
equations.  When  sweeping  in  the  opposite  direction,  the  algorithm  remains  stable  for 
two  to  five  iterations  and  then  diverges  rapidly.  The  time  step  also  affects  the  rate  of 
divergence.  Alternating  the  sweep  direction  may  present  the  advantage  that  every  second 
cycle,  the  parabolic  nature  of  the  equations  is  respected  in  regions  of  reversed  flow. 
Hence,  if  convergence  of  the  algorithm  is  slowed  down  because  the  reversed  flow  region 
is  swept  against  the  local  flow  direction,  alternating  the  sweep  direction  may  improve 
overall  convergence.  Furthermore,  alternating  the  direction  of  sweep  for  the  entire  field 
simplifies  the  overall  sweeping  scheme. 
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The  last  sweeping  scheme,  CMPGB1D,  uses  the  pressure  distribution  provided  by 
ARC2D  to  break  the  computational  field  in  two  segments  about  the  leading  edge 
stagnation  point.  Separate  solutions  are  obtained  for  the  lower  and  upper  surface 
segments  with  each  segment  being  swept  in  its  general  downstream  direction.  The  two 
upstream  stations,  the  stagnation  and  adjacent  stations,  are  fixed  to  the  Navier-Stokes 
data.  They  are  not  updated,  but  still  influence  the  solution  since  they  are  the  upstream 
boundary  conditions  in 

It  should  be  noted  that  the  current  approach  treats  the  lower  and  upper  wake  separately. 
This  may  lead  to  a  discontinuity  of  the  flow  variables  at  the  wake  centerline.  This  is  not 
considered  a  problem  for  the  present  study  and  could  be  corrected  in  later  versions  of 
GBL.  It  would  involve  use  of  points  on  opposite  sides  of  the  computational  domain  to 
ensure  continuity  at  the  wake  centerline.  Such  treatment  is  not  included  here. 


Equation  Level 

Solution  of  individual  equations  or  the  update  of  state  variables  and  the  turbulence  model 
are  done  at  the  equation  level.  GBL  also  includes  a  scheme  to  update  the  inverse  forcing 
function  xw  to  match  pressure  at  the  edge.  Recall  from  section  3  that,  for  isentropic 
flows,  pressure  and  velocity  are  related  at  the  edge.  The  current  update  scheme  matches 
ue  in  the  reverse  flow  region.  The  five  parts  of  the  equation  level  are  now  described 
separately. 


Momentum  equation 

The  sweep  control  level  sets  the  range  [MINXI,MAXXI]  and  sweep  direction  RSWEEP 
for  the  solution  of  the  momentum  equation.  The  actual  solution  is  done  in  routine 
USWEES.  For  each  station,  the  following  sequence  of  operations  is  required.  First,  the 
velocity  u  is  saved  to  evaluate  the  time  derivative  u,  and  to  compute  the  residuals.  For 
airfoil  cases,  the  next  operation  is  to  set  the  switch  LSMODE  to  true  ox  false.  Selection 
of  the  mode  for  each  station  is  actually  done  as  part  of  the  inverse  mode  region  update 
which  are  discussed  below. 

The  direct  and  inverse  mode  computations  are  handled  by  separate  blocks  of  routine 
USWEES.  Each  block  evaluates  the  matrix  elements  for  the  tridiagonal  or  modified 
tridiagonal  matrix  solvers,  sets  the  boundary  conditions,  and  solves  the  system  of 
equations.  The  user  has  partial  control  over  the  boundary  conditions.  Over  a  solid 
boundary,  Dirichlet  or  zero  gradient  Neumann  (slope  equal  zero)  boundary  conditions 
may  be  applied  at  the  lower  and  upper  rj  limits  of  the  profile.  In  the  wake,  a  zero 
gradient  Neumann  condition  is  imposed  at  both  boundaries  irrespective  of  the  mode.  The 
user  controls  the  boundary  conditions  through  variables  BCDXL  and  BCDXU  for  the 
direct  mode,  and  through  variables  BCIXL  and  BCIXU  for  the  inverse  mode. 
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The  matrix  solver  returns  the  solution  at  time  n+ 1  in  the  vector  DT.  Under  or  over¬ 
relaxation  is  then  applied  to  the  solution  depending  on  the  factor  WX.  This  factor  is 
added  mostly  for  convenience  and  is  set  to  unity  for  most  cases.  Note  that  by  updating 
the  values  of  u  after  each  station,  this  data  influences  the  computations  at  the  next  station. 


Energy  equation 

Solution  of  the  energy  equation  is  similar  to  that  for  the  momentum  equation.  Routine 
HSWEES  solves  the  equation  on  a  per  station  basis  for  the  range  and  direction  set  at  the 
sweep  control  level.  The  first  operation  is  to  save  the  solution  at  time  n.  The  implicit 
and  forcing  term  coefficients  are  then  computed,  the  boundary  conditions  applied,  and 
the  resulting  tridiagonal  system  of  equations  solved.  The  user  selects  between  Dirichlet 
and  zero  gradient  Neumann  conditions  at  the  lower  and  upper  boundaries  through 
variables  BCDEL  and  BCDEU  respectively.  This  selection  applies  to  all  ^  stations 
including  those  in  the  wake.  The  solution  may  be  under  or  over-relaxed  by  setting  the 
variable  WH  to  a  value  other  than  unity. 


Continuity  equation 

The  continuity  equation  simply  relates  the  components  of  velocity  and  density  to  ensure 
that  new  mass  is  not  created.  Using  the  u  components  of  velocity  and  density  obtained 
from  previous  steps  of  the  solution  cycle,  the  6  components  of  velocity  arc  obtained  with 
an  explicit  integration  process.  The  value  of  u  is  required  at  one  location  to  start  the 
integration.  In  the  present  case,  it  is  known  to  be  zero  on  the  body  surface,  due  to  the  no 
slip  condition,  and  it  is  also  zero  on  the  wake  centerline  by  definition.  The  integration 
process  is  carried  out  in  module  CINTEG  over  the  range  [MINXI,MAXXI],  Direction  of 
the  sweep  has  no  effect. 

Although  the  integration  is  explicit,  the  tridiagonal  implicit  solver  is  used  to  solve  each 
profile.  However,  only  the  lower  and  main  diagonals  are  used,  which  yields  a  bidiagonal 
system  of  equations.  This  means  that  matrix  elements  are  computed  for  the  upper 
boundary  of  the  profile  while  the  Dirichlet  boundary  condition  po  =  0  is  applied  at  the 
lower  boundary.  Solution  of  the  resulting  matrix,  using  the  tridiagonal  solver,  yields 
values  of  po  for  the  profile.  The  values  of  v  are  obtained  by  dividing  the  matrix  solution 
by  the  local  density.  Values  of  the  contravariant  velocity  V  are  also  updated  at  the  same 
time.  Note  that  the  other  contravariant  velocity  component,  U,  does  not  depend  on  o  and 
need  not  be  updated. 


State  variables  &  turbulence  model 

The  state  variables,  p  and  pm,  are  properties  of  the  flow  which  depend  only  on  the  local 
level  of  energy,  i.e.  H,  p,  and  the  velocity  components.  Th'  ir  computation  is  done  on  a 
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point-by-point  basis  within  the  range  [MINXI.MAXXI]  and  without  need  to  sweep  in  a 
particular  direction.  The  perfect  gas  law  is  used  in  routine  UPRHO  to  compute  density. 
Routine  UPRMM  uses  Sutherland's  approximation,  with  reference  temperature  to 
compute  molecular  viscosity. 

The  Baldwin-Lomax  turbulence  model  is  used  in  routine  UPRMT  to  compute  turbulent 
viscosity.  This  process  involves  all  points  in  a  given  profile  because  the  maximum 
velocity  differential  must  be  computed.  Only  stations  in  the  range  [MINXI.MAXXI]  are 
updated. 


Inverse  mode  region  treatment 

Routine  SCHECK  is  an  algorithm  to  detect  regions  of  reversed  flow  on  the  airfoil.  The 
airfoil  surface  is  scanned  from  the  lower  trailing  edge  to  the  upper  trailing  edge,  omitting 
the  forward  stagnation  point  (i.e.  a  zone  cannot  include  the  leading  edge  stagnation 
point).  For  each  scan,  temporary  pointers  are  set,  based  on  the  velocity  gradient  at  the 
wall,  to  mark  the  beginning  and  end  of  up  to  four  zones  of  reversed  flow.  These  are 
normally  found  immediately  after  a  shock  wave  or  near  the  upper  trailing  edge  of  a 
transonic  airfoil  at  positive  angie  of  attack.  The  location  of  the  temporary  pointers  is 
then  compared  to  that  of  permanent  pointers  and  adjusted  only  if  a  zone  expands. 
Allowing  the  reversed  flow  region  to  shrink  results  in  oscillations  and  non-convergence 
of  the  code. 

Having  identified  the  regions  of  separation,  routine  UPTAUW  updates  the  inverse 
forcing  function  xw  if  the  user  has  enabled  the  pressure  matching  option  (LSPMAT  = 
true).  This  algorithm  ujpdates  the  wall  shear  stress  according  to  a  formula  proposed  by 
Van  Dalsem  and  Steger'6' 


=VW+  \Xlw\  (ue~UNDT)  (5-J) 

The  inverse  forcing  function  is  not  updated  after  every  cycle  because  a  change  of  xw 
often  produces  perceptible  changes  of  the  velocity  profile  over  several  iterations,  and 
changing  the  forcing  function  too  soon  leads  to  instability  of  the  algorithm.  The  logic  to 
decide  when  to  update  xw  is  incorporated  in  routine  UPFORF.  It  monitors  the  variation 
of  usdt  from  one  cycle  to  the  next,  and  updates  xw  only  if  the  percentage  change  is  less 
than  a  user  supplied  value  UECHK,  but  not  before  three  iterations  have  been  completed. 
This  last  condition  results  from  experience  gained  with  GBL,  i.e.  that  as  a  rough  rule-of- 
thumb,  at  least  three  to  five  iterations  are  required  before  updating  the  forcing  function. 
The  problem  with  this  method  is  that  the  tolerence  UECHK  must  be  fairly  coarse  during 
the  initial  cycles  to  allow  more  adjustments  of  t*,,  but  gradually  tightened  as 
convergence  is  achieved.  To  partially  correct  for  this,  the  user  stipulates  the  maximum 
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number  of  iterations,  MAXCNT,  beyond  which  the  forcing  function  is  updated.  Once 
the  value  of  u^dt  has  converged  to  ue  within  a  user  supplied  tolerance  UETOL,  say  16%, 
routine  UPTAUW  is  disabled. 


Function  Level 

The  matrix  elements  are  evaluated  at  the  function  level,  the  most  basic  level  of  the 
current  numerical  method.  Recall  that  the  finite  difference  molecule  yields  a  tridiagonal 
system  of  equations  in  rpl  It  results  that  at  each  point,  the  momentum,  energy  or 
continuity  equation  is  conveniently  expressed  as 

Ak9j,k-\  +  +  Cke%l  i  =  Dk  (5-2) 

where  0  is  u  for  the  momentum  equation,  H  for  energy,  or  p\j  for  continuity.  A,  B  and  C 
are  the  coefficient  of  the  implicit  side  while  D  is  the  forcing  term.  For  the  modified 
tridiagonal  system  of  equations  used  to  solve  the  inverse  mode  form  of  momentum,  the 
above  equation  becomes 

Fk&l+2X  +  A&ltU  +  BkGntf  +  Ckej£  I  =  Dk  (5-3) 

Basic  functions  are  used  to  evaluate  the  implicit  coefficients  and  the  forcing  term  at  a 
given  £,  station.  Since  the  formulation  of  each  equation  changes  near  the  boundaries  of 
the  computational  field,  different  basic  functions  are  used  at  these  locations.  This  is  the 
case  for  the  two  4  stations  at  each  end  of  the  flow  field.  Different  forms  of  the  implicit 
coefficient  B  and  forcing  term  D  are  used  for  the  direct  and  inverse  mode  of  the 
momentum  equation.  Furthermore,  the  inverse  mode  functions  for  the  momentum 
equation  are  different  for  body  and  wake  stations. 

Use  of  the  basic  functions  makes  the  algorithm  very  flexible.  Changes  to  the  equations 
simply  require  new  basic  functions. 
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CODE  TESTING 

This  section  presents  selected  results  to  demonstrate  the  ability  of  the  code  to  compute 
boundary  layer  flows.  First,  a  flat  plate  grid  geometry  is  used  to  verify  coding  accuracy 
for  the  direct  and  inverse  mode  algorithms.  These  tests  involve  laminar  flow  only.  The 
implementation  of  the  Baldwin-Lomax  turbulence  model  is  then  verified  by  comparing 
flat  plate  computations  against  experimental  results. 

The  second  part  of  this  section  presents  results  from  simple  numerical  experiments  to 
explore  the  convergence  characteristics  of  the  code.  The  effect  of  the  time-like  variable 
and  sweep  direction  are  discussed.  In  the  third  part,  the  ability  of  the  code  to  solve  the 
equations  for  an  airfoil  under  attached  flow  conditions  is  shown.  In  the  last  part,  airfoil 
results  for  more  challenging  conditions,  with  embedded  supersonic  and  reversed  flow 
regions,  are  presented. 


Code  Verification 

The  implementation  of  the  algorithms  is  verified  by  conducting  numerical  experiments 
and  comparing  the  results  to  analytic  solutions  for  the  direct  mode,  and  to  other 
computational  results  for  the  inverse  mode.  Each  test  problem  exercises  one  mode  only. 
The  flow  is  laminar  and  the  grid  is  kept  as  simple  as  possible,  a  flat  plate  in  this  instance. 
Results  from  this  type  of  problems  are  very  informative  with  respect  to  issues  as  diverse 
as  grid  dependencies,  convergence  characteristics  and  accuracy.  The  test  direct  mode 
test  case  can  also  be  used  to  verify  the  implementation  of  the  turbulence  model  by 
computing  a  turbulent  flat  plate  flow  and  comparing  the  results  to  experimental  data. 
Results  for  the  three  test  problems  are  now  presented  and  discussed  while  the 
convergence  issue  is  addressed  later  in  the  section. 


Direct  mode  test 

The  direct  mode  algorithm  is  verified  against  laminar  Falkner-Skan  solutions  because  a 
wide  range  of  pressure  gradients,  from  separation  to  strongly  accelerated  flows,  may  be 
simulated.  Falkner-Skan  flows  are  also  laminar  which  eliminates  uncertainties  associated 
with  the  turbulence  model  and  its  experimental  roots.  For  the  current  case,  freestream 
conditions  are  selected  to  yield  incompressible  laminar  flow: 


Re„  =  116,500 
=  0.005 
Too  =  288  K 
1=  1.0 
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Figure  12.  Falkner-Skan  results  for  three  flow  types 
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and  are  kept  constant  for  all  Falkner-Skan  computations.  Similarly,  a  single  grid  with  90 
streamwise  stations  and  50  points  in  the  normal  direction  is  used  for  all  cases.  It  is 
generated  internally  using  the  parameters: 

ss  =  0.01 
=  1.00 

As0  =0.0001 

ft  =  0.2 

fu  =  2.5 

which  yield  a  streamwise  expansion  factor  es  =  0.0775  and  a  maximum  normal 
expansion  factor  en  =  0.0866  in  equation  (2-11).  Note  that  the  grid  starts  at  s  =  0.01 
instead  of  s=0  because  equations  (4-12)  and  (4-14)  for  the  variables  £  and  u, 
respectively,  are  singular  at  this  last  point.  This  complicates  the  initialization  of  the 
velocity  profiles  locally. 

Three  cases  are  computed  for  values  of  the  self-similar  parameter  m  =  -0.08,  m  -  0.0  and 
m  =0.5.  These  correspond  to  flow  under  adverse  pressure  gradient,  constant  pressure 
flow  (i.e.  flat  plate  flow),  and  accelerating  flow  conditions,  respectively.  The  theoretical 
solutions  are  used  to  initialize  the  flow  field  variables  in  each  case.  To  verify  that  the 
numerical  solution  really  converges  to  the  theoretical  solution,  the  initial  data  is 
corrupted  by  setting  the  velocity  values  u  and  D  at  all  non-boundary  points  to  98%  of  the 
theoretical  value.  The  boundary  conditions  and  upstream  conditions  (at  the  first  two 
stations)  are  set  to  their  correct  theoretical  values.  During  the  solution  process,  the  time¬ 
like  variable  At  is  set  to  100.  Results  for  the  three  test  cases  are  shown  in  Figure  12. 
Excellent  agreement  is  obtained  between  the  numerical  results  and  theoiy  which 
demonstrates  the  accuracy  of  the  direct  mode  algorithm  for  laminar  flow. 

In  achieving  the  above  results,  careful  selection  of  the  grid  is  necessary.  For  flows  subject 
to  a  pressure  gradient,  it  is  important  to  concentrate  sufficient  streamwise  stations  to 
resolve  the  pressure  terms;  otherwise,  the  accuracy  is  degraded.  This  applies  to  stations 
near  the  leading  edge  where  pressure  gradients  are  strong,  and  near  the  trailing  edge  of 
the  flat  plate,  where  spacing  of  the  streamwise  stations  is  coarse.  Furthermore,  near  the 
trailing  edge,  the  coarseness  of  the  grid  even  results  in  different  values  of  the  pressure 
terms  if  a  thre“-point-backward  finite  differences  is  used  instead  of  the  centered 
difference. 

The  extent  of  the  grid  in  the  normal  direction  is  also  important.  If  the  grid  does  not  cover 
the  full  extent  of  the  boundary  layer,  the  boundary  condition  is  imposed  at  a  location 
where  it  does  not  apply  and  affects  the  shape  of  the  local  velocity  profile.  This  problem 
is  experienced  with  the  flat  plate  grid  generation  scheme  built  into  GBL.  It  is  partially 


UNCLASSIFIED 


s/u-  dlVds  u,/u_  C,[Re] 


42 


UNCLASSIFIED 


Figure  13.  Klineberg  flow  results 
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eliminated  with  the  use  of  the  //  and  fu  multiplication  factors,  but  the  fact  remains  that 
the  growth  of  the  boundary  layer  grid  is  that  for  a  constant  pressure  flow.  GBL  would 
benefit  from  an  improved  grid  generation  scheme,  or  a  truly  adaptive  grid  scheme  which 
adapts  itself  to  the  flow  conditions. 

Inverse  mode  test 

The  implementation  of  the  inverse  mode  algorithm  is  verified  by  computing  the  laminar 
flow  over  a  flat  plate  subjected  to  the  linearly  decelerated  flow  of  Klineberg.  Freestream 
conditions  are  the  same  as  those  used  for  the  Falkner-Skan  test  cases,  but  a  different  flat 
plate  geometry  must  be  used  for  compatibility  with  Klineberg’s  results.  The  grid 
generated  has  200  streamwise  stations  and  60  points  across  the  boundary  layer.  The 
parameters  are: 


=  0.03 
se  =  6.375 
Asq  =  0.03 
fi  --  0-5 

fu  =  2.0 


The  number  of  streamwise  stations  and  the  spacing  are  selected  to  yield  approximately 
constant  spacing  in  the  s  direction,  thus  more  stations  are  present  in  the  reversed  flow 
region.  The  flow  field  is  initialized  with  constant  pressure  profiles  ( m  =  0)  scaled  with 
Klineberg’s  edge  velocity  distribution. 

Figure  13  compares  the  numerical  results  to  those  of  Klineberg.  The  skin  friction 
distribution  used  to  drive  the  solution  is  shown  in  figure  13(a)  while  figures  13(b)  and 
13(c)  show  the  corresponding  edge  velocity  and  pressure  gradient  results.  The  computed 
edge  velocity  matches  the  Klineberg  results  very  closely,  which  is  noteworthy  since  the 
momentum  equation  uses  a  Neumann  boundary  condition  at  the  upper  r\  edge.  The 
pressure  gradient,  defined  as 


~s  due 

m  =  — - — 

ue  ds 

differs  in  the  reversed  flow  region.  However,  comparison  of  the  current  results  with 
those  of  Klineberg  is  difficult  because  the  above  pressure  gradient  derivative  is  computed 
from  the  numerical  results.  In  particular,  the  evaluation  of  the  term  due/ds  depends  on 
the  grid  and  the  grid  used  in  the  current  method  is  very  different  grid  than  the  one  used 
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Figure  14.  Streamlines  for  Klineberg  flow 


Figure  15.  Turbulent  flat  plate  flow  -  Andersen  experiment 
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by  Klineberg;  a  body  conformal  grid  with  clustering  of  the  T)  lines  near  the  wall  is  used 
here,  while  Klineberg  uses  a  uniform  grid  with  coarse  spacing  since  his  method  exploits 
transformed  variables  instead  of  scaled  variables. 

Figure  14  presents  the  streamlines  computed  for  the  current  computations.  They 
correspond  closely  to  the  results  of  Klineberg  and  show  that  GBL  captures  the  salient 
features  of  the  flow. 

Turbulence  model  test 

A  test  of  the  turbulence  model  implementation  is  carried  out  by  computing  the  flat  plate 
flow  measured  by  Andersen^1 5\  The  experimental  conditions  for  this  flow  are: 

Re„  =  6.5x10s 
M  >  =  0.028 
=  293  K 
1=  1.0 

Andersen  used  a  wind  tunnel  with  a  2.44  meter  long  working  section  to  simulate  flat 
plate  flow  and  measured  the  mean  velocity  profiles  at  various  streamwise  stations.  In  the 
current  computations,  the  flat  plate  is  approximated  by  a  grid  with  90  streamwise 
stations,  50  points  across  the  boundary  layer  and  the  following  parameters: 

ss  =  0.01 
se  —  2.44 

Aso  =  0.005 
//  =  0-5 
A  =  1.3 


The  flow  field  is  initialized  with  the  laminar  flow  variables.  The  momentum  equation  is 
solved  with  a  Dirichlet  condition  at  the  upper  boundary  in  rj.  Andersen  presents  his 
results  in  boundary  layer  coordinates  u+  and  n+  which  are  defined  as 


u 


+ 


(6-4) 
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Figure  16.  Dependence  of  convergence  on  the  magnitude  of  the  forcing  function  term 


Figure  17.  Effect  of  the  momentum  equation  solver  and  viscous  terms  on  convergence 
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n 


+ 


uxti 

v 


(6-5) 


Andersen’s  results  are  for  Ree  =  2570,  where  Ree  is  the  Reynolds  number,  based  on  the 
momentum  thickness.  In  figure  15,  the  GBL  results  are  compared  for  two  profiles.  For 
the  first  profile  the  Reynolds  number  is  matched  while  for  the  second  profile,  the  value  of 
u+  is  matched  at  the  upper  edge  of  the  flow.  Fairly  good  agreement  is  obtained  between 
the  experimental  data  and  the  GBL  results. 


A  Word  on  Convergence 

The  convergence  characteristics  of  the  boundary  layer  equations  are  difficult  to  quantify 
due  to  their  non-linear  character,  the  number  of  factors  involved  and  the  coupling 
between  some  of  these  same  factors.  The  author  has  nevertheless  gained  sufficient 
experience  with  GBL  to  identify  the  major  convergence  characteristics  of  the  code.  Two 
basic  mechanisms  are  at  work.  The  first  one  is  a  convection  mechanisin  which  acts  while 
a  particular  equation  is  being  solved  to  update  a  given  variable,  e.g.  u  .  At  any  given 
station,  solution  of  the  variable  is  affected  by  the  upstream  values  of  u  since  the  finite 
difference  to  its  derivative  involves  upstream  values  of  u.  Recall  that  while  solving  this 
equation,  all  other  variables  are  kept  constant;  only  the  values  of  the  u  field  are  modified. 
The  second  mechanism  is  more  difficult  to  identify.  It  involves  the  interplay  between  the 
different  equations,  i.e.  how  the  solution  of  the  momentum  equation  affects  the  solution 
of  the  energy  equation,  and  so  on. 

Both  mechanisms  are  affected  by  the  sign  and/or  magnitude  of  the  convective,  viscous 
and  time-step  terms.  Furthermore,  the  convective  and  viscous  terms  are  coupled  while 
the  effect  of  the  time-step  terms  is  separate.  Examples  are  now  presented  and  discussed 
to  indicate  the  role  of  the  above  factors. 


Convective  and  viscous  terms 

Convective  and  viscous  terms  depend  strongly  on  the  streamwise  momentum  flow  field. 
The  former  are  function  of  the  magnitude  of  the  streamwise  velocity  while  viscous  terms 
are  function  of  the  normal  gradient  of  the  same  quantity.  As  a  result,  the  momentum 
equation  is  the  driving  force  behind  the  overall  iterative  scheme.  The  evolution  of  the 
momentum  flow  field  is  itself  driven  by  the  forcing  function,  whether  in  the  direct  or 
inverse  mode.  Figure  16  shows  the  effect  of  the  forcing  term  on  convergence.  It  is 
observed  that,  at  least  for  Falkner-skan  flows,  accelerated  flows  converge  faster  than 
decelerated  flows.  Separated  flows  are  even  slower  to  converge  as  demonstrated  by  the 
solution  of  the  Klineberg  flow;  400  iterations  reduces  the  residual  by  only  five  orders  of 
magnitude  (the  same  time-step  is  used  for  all  computations).  Based  on  this  result,  one 
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Manifestation  of  the  convergence  mechanism  between  the  various  equations 


Figure  19.  Effect  of  the  direction  of  sweep  on  convergence  -  part  1 
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could  speculate  that,  for  airfoil  flows,  the  regions  of  accelerated  flow  converge  faster  than 
regions  of  decelerated  or  separated  flow.  Hence,  the  overall  convergence  is  governed  by 
the  regions  of  decelerated  flow. 

Recalling  that  different  forcing  functions  are  used  for  attached  and  reversed  flow  regions, 
hence  how  do  we  know  if  the  solver  affects  the  convergence?  Figure  17  answers  this 
question.  The  curves  labelled  direct  and  inverse  show  the  convergence  histories  for  the 
same  flat  plate  flow  solved  using  the  direct  and  inverse  solvers  respectively.  The  two 
curves  have  similar  convergence  rates  which  supports  the  hypothesis  that  it  is  the 
magnitude  and  sign  of  the  forcing  function  which  drives  convergence,  not  the  method 
used  to  solve  the  momentum  equation.  The  third  curve  of  Figure  17  shows  the 
convergence  for  the  computation  of  the  Andersen  flow.  It  uses  a  different  grid  than  the 
other  curves,  but  it  has  the  same  time-step.  It  is  included  to  demonstrate  the  effect  of  the 
turbulent  terms. 

The  slope  of  the  turbulent  curve  is  less  than  for  the  laminar  curve,  which  indicates  that 
the  turbulent  terms  decrease  the  convergence  rate  of  the  algorithm.  These  terms 
contribute  to  the  magnitude  of  the  elements  of  the  left  hand  (or  implicit)  side  of  the 
equation  only.  The  relative  magnitude  of  the  forcing  term,  located  on  the  right  hand  (or 
explicit )  side  of  the  equation,  is  therefore  decreased  and  the  overall  convergence  is 
slower. 

Another  factor  may  be  responsible  for  the  slower  convergence  of  the  turbulent  flow  case 
with  respect  to  the  laminar  flow  cases:  the  quality  of  the  initial  guess.  Assuming  that 
convergence  is  rate  limited,  which  it  appears  to  be  for  a  given  time  step  and  forcing 
function,  more  iterations  are  required  to  decrease  the  residual  to  machine  zero  when  the 
initial  guess  is  further  away  from  the  solution.  For  the  turbulent  case  of  figure  17,  the 
initial  guess  is  a  laminar  solution,  which  is  quite  different  from  the  converged  solution. 

Figure  18  shows  the  typical  relation  between  the  momentum,  energy  and  density 
residuals.  The  three  variables  converge  to  different  values,  which  reflects  the  sensivity  of 
each  variable  on  the  others,  i.e.  energy  and  density  are  less  sensitive  to  changes  of  the 
velocity  than  velocity  to  changes  of  density  or  energy.  More  importantly,  each  equation 
requires  a  different  number  of  iterations  o  converge  with  the  convergence  of  the 
momentum  flow  field  lagging  that  of  the  energy  and  density  fields.  This  is  a  clear 
indication  of  the  presence  of  a  mechanism  oetween  the  various  equations. 

Alternating  the  direction  of  sweep  produces  a  sawtooth  convergence  history  as  illustrated 
in  figure  19.  The  scheme  converges  while  sweeping  in  the  downstream  direction,  but 
diverges  in  the  upstream  direction,  at  least  for  the  time-step  used.  In  fact,  it  is  the 
experience  of  the  author  that  the  scheme  is  unconditionally  stable  when  sweeping  in  the 
downstream  direction,  but  only  conditionally  stable  when  sweeping  in  the  upstream 
direction. 
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Figure  20.  Effect  of  the  direction  of  sweep  on  convergence  -  part  2 
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For  the  flat  plate  case,  the  stability  margin  lies  near  the  unit  time-step  value.  Figure  20 
compares  the  effect  of  downstream  and  upstream  sweeps  for  this  case.  The  downstream 
case  converges  slowly  while  the  upstream  case  shows  an  increasing  rate  of  divergence 
over  the  first  fifty  iterations,  and  then  levels  out  to  a  large  value.  The  flow  field  variables 
display  an  oscillatory  behaviour.  For  smaller  time-steps,  the  scheme  converges,  albeit 
very  slowly  (thousands  of  iterations)  and  at  the  cost  of  a  loss  in  accuracy,  as  will  be 
discussed  in  the  "Time-step  terms"  subsection  below. 

Alternating  the  direction  of  sweep  does  not  improve  the  convergence  of  the  Klineberg 
case  either.  The  reversed  flow  region  is  so  small  with  respect  to  the  attached  flow  field 
that  any  local  convergence  improvements  cannot  influence  the  overall  convergence. 
Hence,  there  is  no  gain  from  alternating  the  direction  of  sweep.  Then,  why  is  the 
reversed  flow  region  converging  at  all  when  all  sweeps  are  in  the  downstream  direction 
and  the  time-step  is  large?  The  answer  resides  in  the  magnitude  of  the  contravariant 
velocity  U.  As  a  general  rule,  the  magnitude  of  U  is  small  in  the  separation  bubble,  thus 
the  influence  of  the  convective  terms  is  reduced  locally.  One  may  speculate  that  it  is  this 
characteristic  which  allows  the  overall  scheme  to  converge. 

Time-step  terms 

The  time  term  splits  in  two  components,  one  on  the  implicit  side  of  the  equation,  and  one 
on  the  explicit  side.  On  the  implicit  side,  it  increases  the  diagonal  dominance  of  the 
system  of  equations.  On  the  exglicit^side,  it  impresses  the  value  of  the  flow  field  variable 
at  the  previous  time  step  (i.e.  u  ,  H  or  (po)n)  on  the  forcing  function.  For  downstream 
sweeps  of  the  laminar  flat  plate  flow,  the  net  effect  of  decreasing  the  time-step  is  to  slow 
convergence,  as  shown  in  figure  21.  As  the  magnitude  of  the  time-step  is  decreased  by 
orders  of  magnitude  from  l(r  to  1(T3,  more  and  more  iterations  are  required  to  converge. 
Ultimately,  a  lower  limit  for  the  time-step  is  reached  beyond  which  the  accuracy  is 
degraded  because  the  time  term  overwhelms  the  forcing  term  and  truncation  error  sets  in. 


Airfoil  Case  -  Attached  Flow 

The  ability  of  the  code  to  compute  attached  flow  conditions  over  an  airfoil  is  examined 
by  computing  the  flow  over  a  NACA  0012  at  =  0.7  and  a=  1.49°.  This  test  case, 
used  by  Holst^1^  to  compare  several  computer  codes,  produces  attached  flow  over  the 
full  extent  of  the  airfoil.  It  is  mildly  transonic  with  a  maximum  local  Mach  number  of 
1.04  near  the  17%  chord  location. 

The  present  results  are  for  the  standard  configuration  of  the  code  with  flow  checks  to 
detect  regions  of  reversed  flow,  but  without  the  algorithm  to  update  the  wall  shear  stress 
distribution  in  regions  of  reversed  flow.  The  boundary  layer  grid  is  formed  from  a  subset 
of  the  ARC2D  grid  and  the  initial  flow  data  is  from  a  converged  ARC2D  solution. 
However,  the  velocity  data  is  multiplied  by  0.98  to  verify  that  GBL  recovers  the  solution. 
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Downstream  sweeps  are  used  to  solve  the  equations  and  the  solution  is  limited  to  the  top 
surface  of  the  airfoil  in  the  interest  of  clarity.  Note  that  the  lower  surface  is  easier  to 
solve  and  does  not  generate  new  knowledge  about  the  code.  Finally,  only  airfoil  stations 
are  used  since  the  r|  =  1  line  does  not  overlap  the  wake  centerline. 

For  the  first  set  of  results,  the  ARC2D  surface  pressure  is  used  as  the  forcing  function  to 
GBL.  The  edge  velocity,  used  as  the  upper  boundary  condition  for  the  momentum 
equation,  is  computed  from  the  velocity-pressure  relation  of  equation  (4-4).  This  case 
converges  smoothly  to  machine  zero  in  125  iterations.  Figure  22  compares  the  resulting 
velocity  profiles  to  the  ARC2D  solutions  at  15,  25,  35,  50,  75  and  95%  of  chord.  Very 
good  agreement  is  obtained  for  all  stations  except  near  the  trailing  edge  where  there  is  a 
small  difference  in  edge  velocity.  This  occurs  because  the  code  detects  flow  reversal 
over  a  very  small  range  near  the  trailing  edge,  and  switches  to  the  inverse  flow  solver  at 
these  stations. 

The  reversed  flow  region  is  predicted  because  an  appreciable  normal  pressure  gradient 
exists  near  the  trailing  edge,  but  it  is  not  accounted  for  by  GBL.  To  correct  for  this 
condition,  the  pressure  is  interpolated  from  ARC2D  at  the  approximate  boundary  layer 
edge.  However,  how  do  we  find  the  edge?  To  answer  this,  we  must  consider  the 
vorticity  field  around  the  airfoil.  Recall  from  the  airfoil  initialization  section  that  vorticity 
is  present  only  within  the  boundary  layer  and  that  its  intensity  is  largest  near  the  surface. 
As  we  move  out  of  the  boundary  layer  in  the  normal  direction,  it  decreases  by  several 
orders  of  magnitude  until  it  vanishes.  Figure  23  shows  this  variation  for  the  present  test 
case.  It  is  observed  that  the  curve  with  the  value  1 0)  I  =  1  corresponds  closely  to  the 
boundary  layer  edge.  The  resulting  pressure  distribution  is  compared  to  the  surface 
pressure  in  figure  24.  The  normal  pressure  variation  near  the  trailing  edge  shows  clearly. 
GBL  does  not  predict  separation  with  the  new  pressure  distribution  and  the  velocity 
profile  results  near  the  trailing  edge  are  improved. 

In  another  test  run,  the  airfoil  flow  is  solved  using  the  inverse  mode  at  all  stations.  The 
wall  shear  stress  forcing  function  is  obtained  from  a  three  point  approximation  of  the 
Navier-Stokes  data  near  the  surface.  This  solution  requires  35  more  iterations  than  the 
direct  mode  case  and  the  resulting  velocity  profiles  are  not  accurate.  Figure  25  shows 
that  the  edge  velocity  distribution  is  not  captured.  When  run  with  the  edge  velocity 
matching  option,  the  same  case  converges  to  the  correct  solution  in  300  iterations.  The 
larger  number  of  iterations  is  required  because  each  adjustment  of  the  forcing  function  is 
equivalent  to  restarting  the  computations;  an  acceptable  wall  shear  stress  distribution  is 
obtained  in  about  150  iterations.  The  matched  edge  velocity  is  within  half  of  a  percent 
and  the  resulting  velocity  profiles  are  essentially  identical  to  those  for  the  direct  mode 
case.  The  final  wall  shear  stress  distribution,  shown  in  figure  26,  differs  only  slightly 
from  the  initial  distribution,  mostly  in  the  supersonic  bubble  and  immediately  aft  of  the 
shock. 
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Figure  24.  Comparison  of  surface  pressure  to  pressure  along  vorticity  line  1 0)1  =  1 
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Figure  25.  Direct  and  inverse  mode  values  of  edge  velocity 


Figure  26.  Initial  and  final  wall  shear  stress  distributions  for  match  of  edge  velocity 
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Figure  27.  Velocity  profiles  for  adaptive  grid,  direct  mode 
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In  the  above  inverse  mode  test  case,  one  must  be  careful  in  selecting  the  pressure 
distribution  to  be  matched.  The  wall  shear  stress  forcing  function  is  derived  by  applying 
the  momentum  equation  at  the  surface  and  therefore,  it  represents  surface  pressure. 
Using  the  pressure  distribution  of  the  vorticity  line  1 0)1  =1,  even  though  it  is  almost 
identical  to  the  surface  pressure  (see  figure  24),  results  in  a  substantial  change  (=25%)  in 
the  computed  wall  shear  stress. 

The  above  direct  mode  cases  were  run  on  a  subset  grid  from  ARC2D,  but  an  "adaptive” 
grid  can  also  be  used.  Such  a  grid  is  generated  by  GBL  to  overlap  the  ARC2D  grid  with 
the  parameters  Anmin  =  10-5,  fi  =  1.0  and  fu  =  3.5.  The  maximum  number  of  points 
across  the  boundary  layer  is  45.  Using  this  grid  with  the  pressure  distribution  on  the 
I  to  I  =  1  line,  GBL  switches  to  the  inverse  mode  for  all  stations  aft  of  60%  chord 
immediately  upon  completing  the  first  iteration.  As  a  result,  the  edge  velocity  is  not 
matched  for  these  stations,  but  the  velocity  profiles  are  nevertheless  acceptable. 

The  early  switch  to  the  inverse  mode  is  traced  to  the  corrupted  data  which  yields  poor 
estimates  of  the  derivatives  initially.  It  is  corrected  with  a  slow  start  procedure,  i.e. 
under-relaxation  factors  are  applied  to  the  equations  over  the  first  ten  iterations.  For  the 
first  iteration,  the  under-relaxation  factor  is  10%,  and  it  is  increased  by  10%  upon  each 
successive  iteration.  Figure  27  presents  the  resulting  velocity  profiles.  They  agree  very 
well  with  those  of  ARC2D,  although  the  normal  pressure  gradient  effect  can  still  be  seen 
near  the  trailing  edge.  These  results  indicate  that  an  improved  adaptive  grid  scheme  may 
benefit  the  code  since  all  field  points  could  be  used;  however,  under-relaxation  is 
required  to  start  the  computations  to  avoid  tripping  the  inverse  mode.  The  broad  range  of 
flow  conditions  may  also  make  it  difficult  to  find  suitable  criteria  to  adapt  the  grid. 


Airfoil  Case  -  Separated  Flow 

Two  separated  flow  test  cases  are  now  presented.  First,  a  NACA  0012  airfoil  at  3° 
angle-of-attack  and  M ^  ~  0.70  is  considered.  This  flow,  previously  computed  by 
Maksymiuk  et  al.[17',  accelerates  rapidly  over  the  leading  edge  to  create  a  supersonic 
bubble  with  a  maximum  Mach  number  M  =  1.2.  It  decelerates  through  a  shock  at  =33% 
chord  and  immediately  separates  at  the  base  of  the  shock.  The  separation,  however,  is 
small  and  the  flow  reattaches  immediately  up  to  a  point  near  the  trailing  edge  where  it 
separates  again.  It  is  not  clear  that  the  trailing  edge  separation  exists  physically;  it  could 
be  an  artifact  due  to  the  local  skewness  of  the  grid.  Finally,  transition  from  laminar  to 
turbulent  flow  is  fixed  at  5%  chord. 

A  subset  of  Maksymiuk’s  Navier-Stokes  grid  is  used  as  the  boundary  layer  grid.  The 
maximum  number  of  points  allowed  through  the  boundary  layer  is  40.  The  computations 
are  restricted  to  the  airfoil  upper  surface  with  all  sweeps  done  in  the  downstream 
direction.  A  Dirichlet  condition  is  used  at  the  upper  T|  boundary  of  the  momentum 
equation  in  direct  mode  regions  while  a  Neumann  condition  is  used  for  the  inverse  mode. 
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Figure  28.  Velocity  profiles  for  Maksymiuk  case,  velocity  matching  not  enabled 
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Velocity  matching  is  not  enabled  in  the  inverse  regions  and  the  forcing  pressure 
distribution  is  from  the  vorticity  level  I  (hLUt  I  =  1.0. 

The  velocity  profiles  for  the  15,  25,  33,  50,  75  and  95%  chord  stations  are  shown  in 
figure  28  while  figure  29  presents  the  corresponding  edge  velocity  results.  The  edge 
velocity  distribution  is  very  close  to  the  direct  mode  distribution,  but  it  is  still  not  fully 
captured  around  the  shock  and  results  in  significant  differences  in  the  local  velocity 
profiles.  GBL  still  predicts  good  velocity  profiles  for  most  stations.  Use  of  velocity 
matching  in  the  inverse  mode  region  results  in  a  decrease  of  the  xw  forcing  function  aft  of 
the  shock  (see  figure  30)  and  a  velocity  defect  of  the  local  profiles  as  shown  in  figure  31. 
However,  GBL  recovers  well  and  velocity  profiles  downstream  of  the  shock  compare 
well  with  those  from  the  Navier-Stokes  solution. 

The  matching  process  is  slow  and  requires  much  user  interaction.  For  this  case,  it  is  even 
doubtful  that  it  fully  works  since  it  was  required  to  relax  the  matching  criteria  UETOL  to 
1%  to  obtain  convergence.  In  other  words,  the  matched  solution  is  not  much  better  than 
the  one  without  velocity  matching.  From  these  results,  it  appears  that  GBL  has  difficulty 
in  handling  shock  waves.  The  code  may  have  much  to  gain  from  the  application  of  a 
Reimann  problem  at  the  shock,  but  this  matter  is  not  addressed  in  this  work. 

To  further  investigate  the  behaviour  of  GBL  near  shock  waves,  a  second  separated  flow 
case  is  used.  Proposed  by  Holst116*,  it  involves  a  NACA  0012  airfoil  at  8.34°  angle-of- 
attack  and  a  freestream  Mach  number  =  0.55,  which  results  in  a  very  strong 
acceleration  of  the  flow  around  the  leading  edge  to  a  maximum  Mach  number  M  -  1.34, 
and  a  strong  shock  at  13%  chord.  The  deceleration  of  the  flow  through  the  shock  induces 
a  substantial  separation  bubble  at  its  base,  followed  by  reattachment  and  a  second 
separation  zone  at  the  trailing  edge.  The  second  separation  bubble  is  due  to  the  strong 
adverse  pressure  gradient  present  on  the  top  surface  of  the  airfoil.  Laminar  to  turbulent 
transition  is  fixed  at  5%  chord  for  the  computations. 

To  compute  this  flow,  GBL  is  used  without  the  velocity  matching  option  and  with  the 
standard  reverse  zone  detection  algorithm.  Dirichlet  and  Neumann  boundary  conditions 
are  imposed  at  the  upper  T|  limit  of  the  momentum  equation  for  the  direct  and  inverse 
modes,  respectively,  and  the  boundary  layer  grid  is  generated  from  a  subset  of  the 
Maksymiuk  Navier-Stokes  grid  with  a  maximum  of  45  points  in  the  normal  direction. 
The  vorticity  level  lcoCH/l  =0.5  is  used  to  determine  the  edge  of  the  boundary  layer. 
This  value  was  selected  after  examination  of  the  vorticity  contours  of  the  ARC2D 
solution,  as  depicted  in  figure  32.  Note  how  the  shock  wave  causes  a  large  increase  of 
the  boundary  layer  thickness  and  a  corresponding  increase  of  the  magnitude  of  vorticity. 

The  upper  edge  velocity  is  set  from  the  vorticity  level  I  c ocui  I  =  0.5,  but  surface  pressure 
is  used  as  the  driving  force.  It  is  interesting  to  see,  in  figure  33,  how  much  pressure 
varies  across  the  boundary  layer  near  tiie  shock.  The  normal  pressure  gradient  is  already 
significant  twenty  points  above  the  surface,  which  is  still  well  within  the  boundary  layer. 
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Maksymiuk  case 
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Figure  32.  Vorticity  contours  for  NACA  0012  airfoil  at  a  =  8.34°, 
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This  pressure  variation  is  responsible  for  the  convergence  problems  of  GBL  near  strong 
shocks  and  special  handling  of  the  shock  region  is  required.  This  conclusion  is  also 
reached  because  GBL  converges  well  and  produces  velocity  profiles  in  good  agreement 
with  the  Navier-Stokes  solution  when  it  is  run  with  only  stations  ahead  of  the  shock,  or 
only  station  aft  of  the  shock. 

When  the  shock  region  is  included,  GBL  produces  the  velocity  profile  results  of  figure 
34.  Note  how  a  negative  edge  velocity  is  obtained  for  profiles  near  the  shock  wave. 
Recall  that  the  surface  pressure  is  used  to  drive  the  solution,  hence  this  is  a  direct  result 
of  the  normal  pressure  variation.  These  major  changes  to  the  profiles  in  the  shock  region 
propagate  downstream.  GBL  recovers  somewhat  from  this  condition,  but  never 
completely;  as  it  now  thinks  that  the  boundary  layer  edge  is  well  within  the  grid  limits. 
Conversely,  it  is  possible  that  the  grid  is  too  large  and  the  edge  velocity  condition  is  then 
taken  well  outside  the  range  of  validity  of  the  boundary  layer  equations.  This 
observation  is  supported  by  the  good  match  of  the  GBL  profiles  with  the  Navier-Stokes 
solution  near  the  wall  for  most  stations  aft  of  the  shock  location.  The  profiles  in  the 
trailing  edge  region  disagree  completely  with  the  Navier-Stokes  results,  but  this  is  not  too 
surprising  since  the  upstream  velocity  defect  is  compounded  with  a  second  separation 
zone  in  that  region. 

Interesting  results  are  obtained  when  GBL  is  run  with  a  grid  made  up  of  only  the  first 
twenty  grid  points  in  the  normal  direction.  The  edge  velocity  at  the  twentieth  point  is 
imposed  as  a  Dirichlet  boundary  condition  on  the  momentum  equation  for  both  the  direct 
and  inverse  mode  regions.  The  surface  pressure  is  still  used  as  the  driving  function. 
Figure  35  shows  good  agreement  with  the  Navier-Stokes  solution,  even  in  the  shock  and 
trailing  edge  regions  which  are  solved  with  the  inverse  mode.  Use  of  GBL  in  this 
manner  will  be  acceptable  when  it  is  integrated  with  the  Navier-Stokes  code  to  produce  a 
Fortified  Navier-Stokes  procedure. 
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Figure  34.  Velocity  profiles  for  separated  flow  case  of  Holst  -  regular  algorithm 
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Figure  35.  Velocity  profiles  for  separated  flow  case  of  Holst  -  algorithm  limited  to 
twenty  points 
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CONCLUSIONS 

The  current  work  is  aimed  at  achieving  two  principal  goals.  First,  to  develop  a  flexible 
research  tool  to  study  the  behaviour  of  the  Steger  and  Van  Dalsem  algorithm  for  solving 
the  boundary  layer  equations,  and  second,  to  investigate  how  the  same  algorithm  may  be 
used  in  conjunction  with  a  Navier-Stokes  procedure  to  study  the  Fortified  Navier-Stokes 
concept. 

With  respect  to  the  first  goal,  GBL  is  now  developed  to  the  point  of  being  a  useful 
research  tool.  This  report  has  outlined  the  main  features  of  the  code  pertaining  to  its 
input  requirements,  its  various  initialization  procedures,  its  solution  techniques  and 
miscellaneous  computational  tools  required  by  the  code.  Verification  results  presented  in 
the  first  part  of  the  CODE  TESTING  section  confirm  that  the  algorithm  is  implemented 
correctly  and  that  it  is  accurate  for  attached  and  separated  laminar  flow.  It  is  also  found 
that,  at  least  for  these  simple  flows,  the  algorithm  converges  faster  when  subjected  to  a 
negative  pressure  gradient  and  that  it  is  slowest  in  regions  of  reversed  flow.  These  results 
also  show  that  the  rate  of  convergence  of  the  algorithm  depends  on  the  magnitude  of  the 
forcing  function,  and  not  on  the  method  used  to  solve  the  equations,  i.e.  standard 
tridiagonal  solver  versus  modified  tridiagonal  solver. 

Simple  numerical  experiments  with  laminar  flows  show  the  algorithm  to  be 
unconditionally  stable  with  respect  to  the  time  step  used  to  relax  the  equations,  provided 
the  flow  field  is  swept  in  the  general  downstream  direction  and  that  regions  of  reversed 
flow  are  weak.  The  algorithm  is  only  conditionally  stable  when  the  flow  field  is  swept  in 
the  opposite  direction.  Small  time  steps  produce  convergence,  but  the  number  of 
iterations  required  to  obtain  a  converged  solution  is  greatly  increased  in  comparison  to  a 
downstream  sweeping  scheme.  It  should  also  be  noted  that  too  small  a  time  step  reduces 
the  influence  of  the  forcing  function  and  degrades  accuracy. 

Comparison  of  GBL  results  against  the  experimental  data  of  Andersen  show  that  the 
Baldwin-Lomax  turbulence  model  is  succesfully  implemented.  It  also  demonstrates  that 
selecting  an  initial  guess  far  from  the  final  answer  affects  convergence  due  to  the 
maximum  rate  at  which  the  method  can  converge  for  a  given  value  of  the  forcing 
function. 

Numerical  results  for  the  computation  of  transonic  airfoil  flows  demonstrate  the  ability  of 
GBL  to  solve  such  flows  accurately  using  the  direct,  inverse  or  mixed  direct/inverse 
modes.  In  particular,  GBL  predicts  accurate  results  for  fully  attached  flow  conditions 
and  for  cases  with  embedded  supersonic  flow  and  a  small  separated  flow  region  as  long 
as  the  shock  wave  is  weak.  Strong  shock  waves  cause  an  appreciable  pressure  gradient 
across  the  boundary  layer  which  is  not  modelled  in  GBL;  hence  the  present  version  of  the 
code  fails  to  compute  such  a  case.  The  normal  pressure  gradient  which  characterizes  the 
trailing  edge  region  may  also  cause  problems  with  GBL,  but  it  is  easily  handled  by 
selecting  pressure  from  a  point  above  the  surface. 
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Significant  progress  has  been  made  towards  achieving  the  second  goal  of  the  present 
work:  to  develop  a  Fortified  Navier-Stokes  code.  The  framework  to  transfer  GBL  results 
to  the  Navier-Stokes  code,  and  vice  versa,  is  now  in  place.  It  is  also  shown  that  GBL 
predicts  the  velocity  profiles  accurately  when  a  Dirichlet  condition  is  applied  at  the  upper 
boundary  of  the  momentum  equation  with  mixed  direct/inverse  modes,  and  the 
computations  are  limited  to  the  portion  of  the  grid  nearest  the  wall  (say  20  points). 

Future  work  will  concentrate  on  the  integration  of  the  GBL  and  Navier-Stokes  codes  and 
the  development  of  suitable  interaction  mechanisms  between  the  two  procedures.  The 
success  of  the  method  also  hinges  on  the  evolution  of  the  boundary  conditions  and 
forcing  functions  in  the  Navier-Stokes  procedure. 
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APPENDIX  A;  Sample  Input  File 


GBL  Inputs 

Value 

Name 

Description 

9000000. 

RE 

Freestream  Reynolds  number.  Re.,  based  on  freestream  velocity  u_ 

1.000 

REFL 

Reference  length  (in  meters)  for  Re. 

0.70 

FSMACH 

Freestream  Mach  number,  Af_ 

255. 

FST 

Freestream  temperature  in  Kelvin 

1 

RSTURB 

Turbulence  model,  0  for  laminar  flow  and  1  for  Baldwin-Lomax  model 

0.05 

TRANSLO 

TRANSLO  and  TRANSUP  indicate  the  chord  stations  where  laminar  to 

0.05 

TRANSUP 

turbulent  transition  if  fixed.  Also  used  for  the  flat  plate  grid  generation 

3 

RSINCO 

Switch  for  initial  conditions:  0  =  FSK,  1  =  Kline  berg,  2,3  =  Navier-Stokes. 

0.10 

OTA 

Variables  DTA  and  DTW  allow  the  selection  of  different  time  steps 

0.10 

DTW 

for  the  airfoil  and  wake  stations. 

200 

URN 

Number  of  iterations  for  run  (1  to  2000) 

0 

RSUEDG 

Type  of  initialization  for  u.  and  p,  respectively,  0  =  from  surface  conditions, 

0 

RSPEDG 

1  =  from  interpolation  at  1  col  level  1  1,2  =  from  rf  line  KBL 

2 

RSRES 

Residue  type,  0  =  none,  1  =  equation,  2  =  lime  difference 

1 

RSGEOM 

Indicate  kind  of  geometry  used.  0  =  flat  plate,  1  =  ARC2D  airfoil  grid 

15 

NXS 

Number  of  streamwise  grid  stations  for  flat  plate  grid  generator 

45 

KBL 

Number  of  normal  points  for  flat  plate  and  adaptive  grid  generators 

0.01 

XSGRD 

First  coordinate  for  flat  plate  grid  generator,  s. 

1.000 

XEGRD 

Last  coordinate  for  flat  plate  grid  generator,  s4 

0.051 

DXSGRD 

Grid  spacing  between  first  and  second  %  stations 

t  00001 

DYMGRD 

Minimum  ri  spacing  for  adaptive  airfoil  grid  in  INIT3 

0.5 

FLGRD 

Multiplication  factor  for  lower  limit  in  grid  generators 

3.5 

FUGRD 

Multiplication  factor  for  upper  limit  in  grid  generators 

1.0 

WCUT 

Cutoff  value  of  1  col  for  interpolation  of  boundary  layer  edge  conditions 

false. 

LSYMIN 

The  value  DYMGRD  is  imposed  on  the  grid  generator  lower  limit 

true. 

— 

LSGEDG 

Used  in  INTT2,  select  to  use  NDT  with  1  tol  criterion  for  fix  grid  at  KBL) 

UNCLASSIFIED 


UNCLASSIFIED 


GBL  Inputs  (Continued) 

Value  Name  Description 

.true. 

LSMODE 

The  direct  mode  is  used 

false. 

LSFIXM 

Only  the  LSMODE  mode  is  used 

false. 

LSPMAT 

GBL  attempts  edge  pressure/velocity  match 

.false. 

LSSLOW 

Slow  start  is  used  over  the  first  10  iterations 

false. 

LSTGFS 

The  stagnation  stations  are  reset  to  Falkner-Skan  profiles 

.false. 

LSPNOR 

The  pressure  varies  in  the  normal  direction 

.true. 

LSCRPT 

The  initial  velocity  is  corrupted  by  CFACT 

0.98 

CFACT 

Multiplication  factor  to  corrupt  initial  velocity  data 

.false. 

LSGRAP 

The  graphic  option  is  enabled 

1 

RSWEEP 

Sweep  direction: -1  =  -  alternating, 1 

.true. 

LSXMOM 

The  x -momentum  equation  is  enabled 

true. 

LSUcUP 

The  contravariant  velocity  field  Uc  is  updated 

true. 

LSENER 

The  energy  equation  (total  enthalpy)  is  enabled. 

true. 

LSURHO 

The  density  is  updated 

true. 

LSCONT 

The  continuity  equation  is  enabled 

true. 

LSURMM 

The  molecular  viscosity  is  updated 

true. 

LSURMT 

The  turbulent  viscosity  is  updated 

0.00 

VM 

Self-similar  parameter  for  Falkner-Skan  test  case 

190 

JMINXI 

Minimum  station  defining  the  computational  field, 

352 

JMAXXI 

Maximum  station  defining  the  computational  field, 

2 

JOFFST 

Offset  about  each  side  of  the  stagnation  point 

1 

JZOFF 

Number  of  stations  added  on  each  side  of  the  inverse  flow  zxxies 

0 

BCDXL 

Boundary  conditions  at  the  lower  and  upper  r,  boundaries  of  the  x-momentum 

equation  in  the  direct  and  inverse  mode,  as  well  as  for  the  energy  equation. 

0  =  Di  rich  let  condition,  1  =  O'*  order  Neumann  condiuon  (slope  is  zero) 

0 

BCDXU 

0 

BCIXL 

1 

BCIXU 

1 

BCDEL 

1 

BCDEU 
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GBL  Inputs  (Concluded) 

Value 

Name 

Description 

1.00 

WX 

Under-  or  over-relaxation  factor  for  the  x-momentum  and 

1.00 

WH 

energy  equations  respectively 

0.0001 

UECHK(4) 

Hfurr  variation  below  which  is  updated  -  for  each  zone 

0.01 

UETOL(4) 

Convergence  tolerance  between  u^dt  and  u,  -  for  each  zone 

10 

MAXUPD(4) 

Maximum  number  of  iterations  between  updates  -  for  each  zone 

6 

JFV 

Number  of  stations  for  output  of  the  velocity  profiles  (top  surface  only) 

0.150 

FVA(O-IO) 

Chord  stations  to  output 
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APPENDIX  B:  GBL  Direct  and  Inverse  Solvers 

GBL  requires  two  matrix  solvers  for  the  direct  and  inverse  mode,  respectively.  The 
direct  mode  solver,  based  on  the  Thomas  algorithm,  solves  a  tridiagonal  system  of 
equations,  i.e.  for  a  matrix  system 


BfL 
ail+ i 


C/l 

&IL+ 1  C/L+l 

ail +2  B/l+2  C/l+2 


*IL 

D,l 

Xf^+i 

dil+\ 

X/L+  2 

D/L+2 

X/u 

B>iu 

where  the  coefficients  A,  B,  C  and  D  are  provided  by  the  calling  routine,  and  X  is  the 
vector  of  unknows  which  is  solved  for.  Upon  returning,  the  solution  X  is  in  the  D  array. 
The  array  of  diagonal  elements  B  is  also  altered.  The  listing  of  subroutine  THOMAS  is 
given  below. 


SUBROUTINE  THOMAS (IL, IU, AA, BB, CC, DD) 


C...  This  routine  is  straight  out  of  the  Anderson  textbook.  It  is  the 

C...  standard  Thomas  algorithm.  Solution  is  stored  in  DD . 

. . . . 

IMPLICIT  REAL *8  (A-H,0-Z) 

DIMENSION  AA ( 1 ) , BB ( 1 ) , CC ( 1 ) , DD ( 1 ) 

C . 

C...  Establish  upper  triangular  matrix. 

LP  =  IL+1 
DO  10  I  =  LP , IU 
R  =  AA(I) /BB ( I - 1 ) 

BB ( I )  =  BB ( I ) -R*CC ( I - 1 ) 

DD (I)  =  DD (I) -R*DD ( I - 1 ) 

10  CONTINUE 

C . 
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C...  Back  substitution. 

DD(IU)  =  DD(IU) /BB(IU) 

DO  20  I  =  LP, IU 
J  =  IU-I+IL 

DD ( J)  =  (DD (J) -CC { J) *DD ( J+l) ) /BB ( J) 

20  CONTINUE 

C I:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::;:::::::  : 
RETURN 
END 

The  inverse  mode  solver,  developed  by  Van  Dalsem  and  Steger,  uses  a  lower  upper 
decomposition  scheme  to  solve  a  system  of  equations  of  the  form 


where  the  added  coefficients  F  is  due  to  the  inverse  forcing  function.  Once  more,  the 
solution  is  returned  in  array  D.  The  Van  Dalsem  and  Steger  algorithm  is  implemented  in 
subroutine  INVTOM  which  is  listed  below. 

c 
c 

c 
c 
c 
c 
c 

c  . 

IMPLICIT  REAL *8  (A-H,0-Z) 

PARAMETER  (KMAX=100) 

DIMENSION 


SUBROUTINE  INVTOM (IL, IU, A, B, C, D, F) 


This  routine  implements  the  solution  algorithm  of  Van  Dalsem  for 
quasi-tridiagonal  matrices.  It  uses  an  upper-lower  decomposition 
followed  by  backsubstitutions  for  the  solution. 
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>  A(KMAX),  B(KMAX),  C (KMAX) ,  D(KMAX), 

>  F (KMAX) ,  G(KMAX),  P (KMAX) ,  R (KMAX) 

C . 

C...  Compute  the  elements  of  the  lower  and  upper  triangular  matrices. 
P(IU)  =  B ( IU) 

R(IU)  =  F ( IU) 

DO  I  =  IU-1 , IL+1 , -1 
G ( I )  =  C(I) /P(I+1) 

P  (I)  =  B  (I)  -G  ( I )  *A  ( 1  +  1 ) 

R(I)  =  F (I) -G (I) *R ( 1+1 ) 

ENDDO 

G  (IL)  =  C  ( IL)  /P  (IL+1) 

P(IL)  =  B (IL) -G (IL) *R ( IL+1) 

C . 

C. . .  Interim  solution  for  multiplication  by  upper  triangular  matrix. 
C. . .  Results  are  stored  in  array  G. 

C(IU)  =  D ( IU) 

DO  I  =  IU-1 , IL, -1 

C(I)  =  D ( I ) -G ( I ) *C(I+1) 

ENDDO 

C . 

C...  Compute  final  solution.  The  results  are  stored  in  array  D. 

D ( IL)  =  C ( IL) /P ( IL) 

D ( IL+1 )  =  (C ( IL+1 ) -R ( IL+1 )*D(IL))/P( IL+1 ) 

DO  I  =  IL+2 , IU 

D ( I )  =  (C(I) -R(I) *D(IL) —A ( I )  *D  ( I  —  1 )  ) /P  (I) 

ENDDO 

C : : : : : : : : : : : :::::: : : : : : :  :  :  :  :  :  :  :  :  :  : : : : :  :  :  :  :  :  :  :  :  :  :  :  ::::::::  :  :  ::::::  :  :  :  :  : 
RETURN 
END 
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